a” ‘
r‘ t
3% "

p"

e...
. Mud

Hr: ’
'.

L
m!

J'
,9;-
.I'
“.1.
1'”,

 

 

{b.5135

MICHIGAN STA‘I'E Ur RSITY LIBRARIES

minim!” ”i“Minimum

31 93 00904 6271

II

 

 

 

 

 

 

 

 

 

 

 

 

 

 

This is to certify that the

dissertation entitled

Stucture and Dynamics of
Low-Dimensional Systems

presented by

Weiqing Zhong

has been accepted towards fulfillment
of the requirements for

Ph. D . degree in Physics

 

 

W/ [
Major professor

e aqn (/31 (A? (993 Professor David Tomanek

 

Dat

 

MSU 13 an Affirmative AcrinnVEquul Opportunity Institution 0‘ 12771

 

LIBRARY
Michigan State
University

 

 

 

PLACE lN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or betore date due.

DATE DUE DATE DUE DATE DUE

__ Ii
I:JL___i
E]

. ll II I

 

 

 

 

 

 

 

 

i

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

—
L

 

 

 

hi i

MSU Is An Aflirmative Action/Equal Opportunity Institution
emana-DJ

 

 

 

STRUCTURE AND DYNAMICS OF
LOW-DIMENSIONAL SYSTEMS

By

WEIQING ZHON G

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the Degree of

DOCTOR OF PHILOSOPHY

Department of Physics and Astronomy

1993

ABSTRACT

STRUCTURE AND DYNAMICS
OF LOW-DIMENSION AL SYSTEMS

by

Weiqing Zhong

I studied structural and dynamical properties of low dimensional systems using
numerical approaches at different levels of sophistication. I focused my interest on
the physical phenomena associated with the interaction of atoms with surfaces. I
studied two topics which are intimately connected with this interaction. One of them
is a procedure to quantitatively interpret atomic force microscopy (AFM) images of
solid surfaces. The other topic is the effect of adsorbate on structural, electronic and

dynamical properties of the substrate.

When studying theory for the AFM, I considered a model system consisting of
a Pd metal tip probing the graphite surface. Using ab initio density functional the—
ory (DFT), I calculated the interaction between this Pd AFM tip and graphite and _.
mapped it onto a parameterized energy functional. The calculation of AFM tip tra-
jectories at different loads revealed that atomic resolution is only achievable for loads
exceeding 5 x 10'9N. Larger loads, on the other hand, are likely to lead to surface
destruction. Calculated variation of potential energy along a trajectory has been used
subsequently to determine atomic scale friction. Assuming no wear and a complete
dissipation of the energy gained by the Pd tip atoms, I found a very small friction

coefficient p z 10" for loads near 10"8 N and an increase in p with increasing load

in agreement with experiments.

I selected H/Pd as the model system to study the effect of adsorbate on the
structural and dynamical properties of the substrate. I based the description of the H-
Pd system on density functional calculations, which were subsequently mapped onto
a parametrized many—body alloy Hamiltonian. Using this Hamiltonian, I calculated
the equilibrium structure of the clean and hydrogen covered Pd (001) and Pd (110)
surfaces, as well as the corresponding surface phonon spectra. The most pronounced
effect of hydrogen is a strong softening of the Rayleigh wave on Pd(001), which is
indirectly related to “hydrogen embrittlement” observed in the bulk. I addressed this
latter problem using molecular dynamics. I studied the equilibrium structure, elastic
properties, and in particular the mechanical breakdown of bulk Pd under tensile
stress, as a function of temperature and hydrogen concentration. My results indicate
that the microscopic origin of “hydrogen embrittlement” is an increased ductility
and plasticity in regions saturated by hydrogen, in agreement with the postulated

Hydrogen Enhanced Local Plasticity mechanism.

ACKNOWLEDGMENTS

I would like to express deep gratitude to my advisor Prof. David Tomanek for
his guidance, understanding, and encouragement during my five years of graduate
study. I thank him for invaluable training in all aspects of research, including research
methods, research manners, organization of ,work, and presentation of results. I will
benefit from this training for the rest of my life. His inspiration and insight highly

motivated my work and made this research enjoyable.

I am grateful to Professor S. D. Mahanti, Professor P. Danielewicz, Professor
R. Tobin, Professor W. Repko, Professor G. F. Bertsch, and Professor W. Pratt for
their service on my Ph.D. guidance committee. I would express my special thanks
to Professor S. D. Mahanti for his careful guidance and advice. I thank Professor

J. S. Kovacs for valuable advice and J. King and S. Holland for constant help.

I would like to thank Professor A. Bulgac, Professor G.F. Bertsch, Professor T. Ka-
plan, Professor P. Duxbury, and Professor M. Thorpe for valuable suggestions and
stimulating discussions. My thanks also go to my dear friends Dr. Y. Cai, Q. Yang,
Dr. Y. S. Li, Dr. G. Overney, Dr. Z. Sun, N. J. Ju, Dr. Y. Wang, J. M. Song,
B. Zhang, B. Lian, Dr. R. Y. Fan, N. Mousseau, J. D. Chen, X. H. Yu, and Q. F. Zhu

for their friendship and help.

I am grateful to my parents and brothers for their love, encouragement, and

support during all my life. My wife Hongjun Yan deserves my deepest gratitude.

I would also like to acknowledge financial support from Michigan State University,

the Office of Naval Research, and the Center for Fundamental Materials Research.

iii

PUBLICATIONS

W. Zhong and D. Tomanek, First-Principles Theory of Atomic-Scale Friction,
Phys. Rev. Lett. 64, 3054 (1990).

W. Zhong, G. Overney, and D. Tomanek, Theory of Atomic Force Microscopy on
Elastic Surfaces, The Structure of Surfaces III, edited by S.Y. Tong, M.A. Van
Hove, X. Xide and K. Takayanagi, Springer-Verlag, Berlin (1991), p243.

W. Zhong, G. Overney, and D. Tomanek, Limits of Resolution in Atomic Force

Microscopy Images of Graphite, Europhys. Lett. 15, 49 (1991).

G. Overney, W. Zhong, and D. Tomanek, Theory of Elastic Tip—Surface Inter-

actions in Atomic Force Microscopy, J. Vac. Sci. Technol. B9, 479 ( 1991).

D. Tomanek, W. Zhong, and H. Thomas, Calculation of an Atomically Mod-
ulated Friction Force in Atomic Force Microscopy, Europhys. Lett. 15, 887
(1991).

D. Tomanek and W. Zhong, Palladium-Graphite Interaction Potentials Based

on First-Principles Calculations, Phys. Rev. B43, 12623 (1991).

W. Zhong, Y. S. Li, and D. Tomanek, E'fl'ect of Adsorbates on Surface Phonon
Modes: H on Pd(001) and Pd(110), Phys. Rev. B44, 13053 (1991).

G. Overney, D. Tomanek, W. Zhong, Z. Sun, H. Miyazaki, S. D. Mahanti, and
H.-J. Giintherodt, Theory for the Atomic Force Microscopy of Layered Elastic
Surfaces, J. Phys. : Cond. Mat. 4, 4233 (1992).

W. Zhong, Y. Cai, and D. Tomanek, Mechanical Stability of Pd—H Systems: A
Molecular Dynamics Study, Phys. Rev. B 46 8099 (1992).

iv

10.

11.

12.

13.

14.

15.

W. Zhong, G. Overney, and D. Tomanek, Structural Properties of Fe Crystals,
Phys. Rev. B 47 95 (1993).

G. Overney, W. Zhong, and D. Tomanek, Structural Rigidity and Low Frequency
Vibrational Modes of Long Carbon Tubules, Z. Phys. D 27 93 (1993).

W. Zhong, D. Tomanek, and George F. Bertsch Total Energy Calculation for
Extremely Large Clusters: The Recursive Approach, Solid State Comm. 86 607
(1993).

W. Zhong, Y. Cai, and D. Tomanek, Computer Simulation of Hydrogen Em-
brittlement in Metals, Nature 362 435 (1993).

G. Overney, W. Zhong, and D. Tomanek, Structural Optimization of Stage—I
Alkali-Metal-Graphite Intercalation Compounds Using Density Functional The-

ory, in preparation.

D. Tomanek, W. Zhong and E. Krastev Stability of Multi-Shell Fullerenes, sub-

mitted for publication.

Contents

LIST OF TABLES iv
LIST OF FIGURES v
1 Introduction 1
1.1 Theory for the Atomic Force Microscopy ................ 2

1.2 Structure and dynamics of H—Pd systems ................ 7
1.2.1 Effect of hydrogen on surface phonon modes of Pd ....... 7

1.2.2 Mechanical stability of Pd-H systems .............. 9

1.3 Structure of the Thesis .......................... 11

2 Theory 13
2.1 Ab initio Density Functional Formalism ................. 13
2.1.1 Density Functional Theory .................... 13

2.1.2 Local Density Approximation .................. 15

2.1.3 Computational details ...................... 17

2.2 Many-Body Alloy Hamiltonian ...................... 20

2.3 Lattice dynamics ............................. 23

vi

2.4 Molecular dynamics . .

3 Theory for the Atomic Force Microscopy

3.1 Calculation of the Pd-graphite interaction ...............

3.2 Limits of resolution in AFM images of gra—phite ............

3.3 First-principles theory of atomic-scale friction .............

3.4 Ideal friction machines

3.5 Conclusions ......

OOOOOOOOOOOOOOOOOOOOOOOOOO

4 Structure and Dynamics for H Loaded Pd

4.1 Application of the Many-Body Alloy Hamiltonian to the Pd—H system

4.1.1 Construction of the Many-Body Alloy Hamiltonian ......

4.1.2 Structure and stability of clean and hydrogen covered Pd(001)

and Pd(110) surfaces .......................

4.2 Phonon structure in the bulk and at the surface of Pd—H systems

i 4.3 Mechanical stability of Pd-H systems ..................

4.3.1 Thermal expansion and the melting transition in Pd ......

4.3.2 Mechanical stability of H-free Pd under tensile stress .....

4.3.3 Mechanical stability of H loaded Pd under tensile stress . . . .

4.3.4 Discussion . . .

4.4 Summary and conclusions ........................

LIST OF REFERENCES

vii

24

29
31
38
48
55.

63

67
68

69

72
80
90
90
96
101
107

111

114

 

List of Tables

4.1 Interaction parameters used in the Many-Body Alloy Hamiltonian for

the Pd-H system .............................. 69

4.2 Relaxations at clean and hydrogen covered Pd surfaces ......... 75

viii

List of Figures

1.1

3.1

3.2

3.3

3.4

3.5

3.6

3.7

3.8

3.9

Schematic picture of the Atomic Force Microscope (AFM) ...... 4

Adsorption energy of Pd on graphite as a function of the adsorption

height ..... ‘ .............................. 32

Dependence of the Pd—graphite interaction energy on the charge den-

sity at the adsorption site ........................ 35
Radial plot of the charge density of a carbon atom ........... 37

Calcualted adsorption energy of Pd on graphite as a function of the

adsorption height ............................. 39
Force between a Pd AF M tip and a graphite surface .......... 42

Atomic-scale corrugation of the monatomic Pd AFM tip scanning graphite

surface ................................... 43

Total charge density of an AF M tip of Pd interacting with the elastic

graphite surface ................... ' ........... 46

Potential energy and force acting on an AF M tip of Pd along a trajec-

tory on the graphite surface ....................... 49

Microscopic friction coefficient as a function of load .......... 52

3.10 Macroscopic friction force as a function of load applied to a large object 54

ix

3.11 Two models for the Friction Force Microscope .............

3.12 Potential energy, friction force, and averaged friction force in the “max-

imum friction microscope” ........................
3.13 Microscopic friction mechanism in the “realistic friction microscope” .

3.14 Average friction force between a Pd tip and graphite, as a function of

the load on the tip and the force constant of the cantilever ......

4.1 Cohesive energy in bulk Pd and PdH as a function of lattice constant,

based on the MBA Hamiltonian .....................

4.2 Surface energy changes for clean and H-covered Pd (001) and Pd (110)

surfaces as a function of first interlayer spacing .............

4.3 Adsorption energy of H on the Pd (001) surface as-a function of the

adsorption height .............................

4.4 Adsorption energy of H on the Pd (110) surface as a function of the

adsorption height .............................

4.5 Phonon dispersion relations for bulk Pd and PdH ...........

57

59

60

64

70

74

77

78

81

4.6 Effect of adsorbed hydrogen on the surface phonon modes of the Pd(001) 84

4.7 Effect of adsorbed hydrogen on the surface phonon modes of the Pd(110) 85

4.8 Equilibrium Pd atomic volume in bulk Pd as a function of temperature,

and same in bulk PdH as a function of hydrogen concentration . . . .

4.9 Atomic trajectories and the pair correlation function in bulk Pd at

different temperatures ..........................

4.10 Distribution of atomic binding energies in bulk Pd ...........

91

95

4.11 Dependence of the deformation of Pd and PdH, on tensile stress . . . 97

4.12 Dependence of the bulk moduli of Pd and PdH; on temperature and

hydrogen concentration a: ......................... 100

4.13 Distribution of binding energies of Pd atoms in bulk Pdes at T =

300 K at zero and critical tensile stress ................. 104

4.14 Distribution of atomic volumes in bulk PdHogs at T = 300 K. . . . . 105

xi

Chapter 1

Introduction

The major task of theoretical physics is to understand the laws which govern na-
ture, use them to explain observed phenomena, and predict the properties of new
physical systems. The subject of interest for theoretical condensed matter physics
are electronic, magnetic, structural, and dynamical properties of solids and liquids.
The most popular technique to address these properties from first principles is the
Density functional theory [Lun 83] in the Local Density Approximation (LDA), which
has been shown to describe ground state properties of solids with a high accuracy.
Since this theory is free of adjustable parameters, I will use it as a basis for all the

calculations in this Thesis.

The large computational requirement associated with ab initio methods such as
the LDA limits the range of its applicability to relatively small systems with a high
symmetry. In order to describe the dynamical behavior of large systems with low
symmetry, I calculate the interatomic forces using a parametrized Many-Body Alloy
Hamiltonian [Zho 91a], which I developed as an efficient interpolation scheme between
ab initio results. In spite of its simplicity, this Hamiltonian describes the interatomic
interactions in simple metals and late transition metals with sufficient accuracy. The

total energy calculations in this Thesis are performed using the combination of the

LDA and MBA formalisms. The dynamical behavior is then calculated using the in-
teratomic forces in the equations of motion. For a canonical ensemble, this procedure
yields lattice dynamics at zero and finite temperatures, and the elastic response to

externally applied loads.

In this thesis, I have applied these calculation tools to study structural and dy-
namical properties of low dimensional systems. I focused my interest on the physical
phenomena associated with the interaction of atoms with surfaces. One important
application of this interaction is the quantitative interpretation of atomic force mi-
croscopy (AF M) images of solid surfaces. For the model system consisting of a Pd
metal tip and graphite. I calculated the AFM images, the limits of nondestructive
imaging, and the origin of atomic-scale friction. The second important topic is the
effect of adsorption on the properties of the substrate. I considered Pd as a model sys-
tem, and studied the effect of adsorbed hydrogen on the surface structure and phonon
spectra. I found these results to be of significance when interpreting bulk phenom-
ena associated with adsorbed hydrogen, such as the “hydrogen embrittlement. In
the following Section 1.1 and 1.2, I will introduce these two topics, emphasizing the

open questions and the general importance of these calculations. The structure Of my

Thesis will be addressed in Section 1.3.

1.1 Theory for the Atomic Force Microscopy

Following its invention by Binnig et al in 1986, the Atomic Force Microscope (AF M)
[Bin 86, Bin 87] has been rapidly evolving into a powerful tool to examine the mor-
phology and local rigidity of conducting and insulating surfaces alike [Tom 89]. The
AFM uses an “atomically” sharp tip to scan the sample surface at a sample-to-tip

separation of a few angstroms. Unlike the more established Scanning Tunneling Mi-

 

croscope (STM) [Bin 82], which is sensitive to the electronic density of states near
the Fermi energy [Ter 83], the AF M probes the force field Fm between an “atom-
ically sharp” tip and the substrate (Figure 1.1). This force Fm is detected by the
deflection of a soft cantilever which supports the tip. During a horizontal scan of the
surface, the AF M measures and records the equilibrium vertical tip position, while
Fe,“ is been kept constant by regulating the deflection of the cantilever. The reso-
lution of atomic-scale features at the surface is possible if the equilibrium tip height
2 is detectably different at inequivalent (e.g. on-top, hollow) surface sites. This is
the case if the corresponding corrugation £322,005 A, which is the sensitivity of the

AFM under ideal conditiOns.

Since the AFM provides real-space image, same as the STM, the ambiguity and
complexity of interpreting diffraction patterns is avoided. The power of the AFM lies
in its ability to detect local structure rather than periodic patterns. These can be
point defects, steps, grain boundaries, and similar on conducting or non-conducting
surfaces. The AF M has been used suCCessfully to image a variety of systems with
atomic resolution. Very good results have been achieved in ionic system such as NaCl
[Mey 90], and layered materials such as graphite [Bin 87]. While the experiments
have primarily focused on achieving atomic resolution, two questions have remained
unanswered so far. The first is, under which conditions atomic-scale features can be
observed by a microscopic tip, and how these features relate to the surface topography.
The second open question, which will be addressed in this Thesis, relates to the

conditions, under which nondestructive imaging can be achieved.

The major challenge for the theory is to determine whether it is possible to achieve
atomic resolution using the AF M on a specific material, under well controlled and

idealized conditions. The potential for atomic resolution is limited in two ways.

 

E. Force
00
Q

Will)
0000
00000

First, if the load applied to the tip, Fe", is very small, the tip-substrate separation
is large, and the corrugation A2 is too small to be observed. For “dull” tips with
multiple apex atoms, the corrugation is further decreased due to the tip/substrate
incommensurability. At large values of Fe“, where the tip is capable of probing
the corrugation of the charge density, the tip-surface distance is very small. Under
these conditions, a destruction of the surface or the tip is likely. The objective of
my calculations is to determine if there is an optimum force, which yields detectable

corrugations on the atomic scale, but which is still non-destructive to the surface.

The AFM has been used to investigate the friction on the atomic scale [Mat 87].
This experiment investigated the friction force experienced by the tip during a scan
with no wear across a perfect surface. In general, friction forces are observed in
nature whenever two bodies in contact are in relative motion [Lan 60]. The force F f

is related to the applied load Fm between the two bodies by
Ff = I‘chr (11)

The friction coefficient a has been found to range typically between 10‘2 for very
smooth interfaces and 1 for rough interfaces. Friction forces have dissipative character:
they couple the macroscopic mechanical degrees of freedom to microscopic degrees
of freedom. In this way, macroscopic mechanical energy is dissipated. At rough
interfaces, dislocation motion and plastic deformations are mainly responsible for the
dissipation of macroscopic energy. This fact makes the friction process with wear
one of the most complex and least understood process in nature. Friction without
wear, on the other hand, occurs between perfect and weakly interacting surfaces and
is much easier to understand. Since in this case, no atomic rearrangement occurs,
it is phonons and electronic excitations which serve as the only source for energy

dissipation. Atomic-scale friction, which can be quantitatively studied using the AFM

operating in a nondestructive manner, falls in this second category. [Sla 93]

The idealized conditions for this atomic-scale friction open the problem to an inde-
pendent fist-principles theoretical study [Zho 90]. In absence of plastic deformations,
friction originates in an atomic-scale corrugation of the interaction potential which is
probed by the tip when scanning the surface under an external force Fm. The two
contributions to this potential are the variation of tip-substrate bond strength and
the work against Fm if the tip-substrate distance varies along the trajectory. The
maximum friction force can be estimated using the variations of the total potential

energy of the tip-substrate system during the scan.

A microscopic description of friction without wear must address the fact that
friction is a non-conservative process. In other words, the friction force depends on
the direction of motion between the two bodies in contact. A closed-loop integral over
such a force yields a nonzero value which corresponds to the dissipated energy. In
other words, the friction force cannot be obtained from a gradient of a potential. A
quantitative study of friction requires the investigation of the microscopic mechanisms
for energy dissipation and of the effectiveness with which macroscopic degrees of

freedom are coupled to microscopic degrees of freedom.

In order to study these process quantitatively, a Friction Force Microscope (FFM)
has been constructed and used to measure the atomic-scale friction [Ove 92]. A very
precise FFM can determine the horizontal and vertical force on the tip simultane-
ously Since different portions of otherwise flat surface may have different friction
coefficients, the FFM may provide the most detailed information and enhanced con-
trast as compared to the AFM image of such a surface. Hence, friction imaging may

prove to be a very useful application of the FFM [Ove 92].

 

7

1.2 Structure and dynamics of H—Pd systems

The interaction of hydrogen with transition metals is a fundamentally interesting topic
with wide ranging technological applications [Ale 78]. The dissociative adsorption of
molecular hydrogen at a metal surface and the subsequent surface and bulk diffusion
of atomic hydrogen is a prototypical surface process. This is a well-suited model
system for the study of the hydrogen-metal bond at different stages of the reaction.
At least of equal importance is the microscopic understanding of the effect of hydrogen
on the structural and dynamical properties of the substrate or host metal. This effect,
which has received less attention in the literature, manifests itself on the surface as
H-induced surface relaxation, surface reconstruction, as well as a change of the surface
phonon spectra. Hydrogen atoms are known to diffuse into many metals easily. In
the bulk, hydrogen can change the equilibrium structure and elastic properties of the
host metal. Technologically, hydrogen atoms in the bulk of transition metals can
also reduce the mechanical stability of the system significantly, which is known as
hydrogen embrittlement. The particular interest in Pd is motivated to a large degree
by the ability of this metal to form hydrides and thereby to act as a medium for

hydrogen storage.

In this Thesis, I will concentrate on studying the effect of adsorbed hydrogen on the
surface phonon modes of Pd, and the microscopic origin of hydrogen embrittlement

in Pd. I will show that both phenomena are intimately related.

1.2.1 Effect of hydrogen on surface phonon modes of Pd

Extensive studies of interaction between adsorbed H and the Pd surface include the
characterization of the H-metal bond, hydrogen induced surface relaxation and re-

construction, and the effect of hydrogen on the the electronic structure of the metal

 

substrate [Tom 86b, Sun 89, Tom 91a, Cat 83, He 88, Beh 80, Rie 83, Rie 84, Tar 86,
N yb 83, Bes 87]. So far, most research effort has been focussed on the adsorption ge-
ometry and on the electronic and structural properties of the adsorption system.
Surface phonons have received far less attention in the recent literature, in spite of
the wealth of information they contain about the nature of bonding at surfaces. This
is caused mainly by the difficulty to measure and calculate reliable surface phonon

dispersion curves throughout the whole surface Brillouin zone.

Only recently, He time-of-fiight spectroscopy [Toe 87, Lab 87, Doa 83, Har 85]
and electron-energy-loss spectroscopy(EELS)[Leh 83, Roc 84, Wut 86, Iba 87, Yos 88]
have been used to measure the dispersion curves of surface phonons on a variety of
metal substrates. The quantitative interpretation of these data is lacking in most
cases, since predictive ab initio calculations (such as “frozen-phonon” calculations)
are computationally very involved. Only in selected cases, Local Density Approxima-
tion (LDA) [Koh 65] calculations have been performed for the high-symmetry modes
[Ho 86].

The majority of published phonon calculations use lattice dynamics based on sim—
ple two- and three-body potentials [Joh 72, Bor 84, Bor 85, Hal 88]. These types of
calculations use bulk and surface interatomic force constants and distances as inde-
pendent parameters which are chosen to fit experimental results [Lah 87]. In general,
these calculations show a good agreement with the observed data. The predictive
power is limited by the generally large number of force constant parameters which
depend on the model, the system and the surface studied. While these calculations
can provide a rough guidance in the interpretation of experimental results, direct

comparisons between different models are of limited use.

More recently, the Embedded—Atom -Method (EAM) [Daw 84, Foi 85, Foi 87,

 

Daw 89, Foi 89] has been used to calculate phonon dispersion relations on surfaces
such as Cu(100), Cu(lll) and Ag(111) [Nel 77, Nel 89, Luo 88]. The major differ-
ence with respect to the calculations quoted above is that all parameters have been
obtained by fitting the measured bulk properties. The EAM has proven to be quite
successful in the prediction of surface phonon spectra and the corresponding changes
of interatomic force constants and distances at surfaces. The major weakness of these
calculations is the limitation to single-component systems, since charge transfers be-
tween different sites are assumed to be zero. Of less importance is the fact that the
success of the EAM technique depends on the type and quality experimental data

and that the fit is based mainly on such observed bulk properties.

To avoid these weak points, I developed a model Many-Body Alloy (MBA) Hamil-
tonian for the Pd-H system, based on ab initio calculations of bulk materials (both
single component systems and alloys). This Hamiltonian can be easily applied to
many alloy systems, and the parametrization can be uniquely determined from a set
of ab initio calculations. This mapping does not leave any adjustable parameters.
The MBA technique is a very useful tool to determine the total energy and forces
acting on individual atoms at the surface and in the bulk. In particular, using lattice
dynamics, the surface phonon modes of clean and hydrogen adsorbed Pd surface can
be calculated. This allows for a detailed analysis of the effect of hydrogen on the

surface phonon spectra of Pd.

1.2.2 Mechanical stability of Pd-H systems

Fundamental understanding of the mechanical stability of transition metals under var-
ious conditions is of great technological importance. Very little is known about the

atomic-level response of bulk metals to applied tensile stress, at varying temperature

10

and hydrogen concentration. The study of these effects is expected to answer some key
questions related to macroscopic materials properties, especially the so-called “hydro-
gen embrittlement”. This hydrogen-induced reduction of the mechanical strength of
metals causes great concerns when considering hydrogen storage in metals, stability

of fusion reactors and underwater structures, and space technology.

In spite of a significant effort to resolve this important problem, there is still sub-
stantial controversy regarding the microscopic origin of “hydrogen embrittlement” in
a given system [Bir 79]. One of the oldest and most commonly referred to mechanisms
of hydrogen embrittlement is the “decohesion mechanism” which associates this effect
with a decreased metal bond strength in the presence of hydrogen [Ste 60, Foi 86]..
The “hydrogen related phase change mechanism” has been suggested as the origin
of “hydrogen embrittlement” in systems where a brittle hydride phase is stabilized
by the presence of hydrogen and the crack tip stress field [Wes 69]. The “Hydrogen
Enhanced Local Plasticity” (HELP) mechanism postulates a hydrogen-induced local

plasticity enhancement at crack tips which facilitates fracture formation [Bea 72].

The most straight-forward way to address the above questions is to use a molec-
ular dynamics (MD) technique. Such a calculation describes the evolution of the
system with time, based on a direct numerical solution of the equations of motion
for individual atoms. This procedure goes beyond lattice dynamics [Mar 63] which
is limited to small atomic displacements and can not address problems such as the
melting transition or fracture. The appropriate MD technique for these questions will
describe the evolution of the system at constant nonzero temperatures using N osé dy-
namics. It can also model the response of a system to an externally applied tensile
stress [All 90, Car 90]. Such a simulation provides microscopic information about the

dynamics of the system, including structural changes during phase transitions (such

11

as melting or crack formation due to tensile stress under different conditions).

1.3 Structure of the Thesis

In this Thesis, I will address primarily two subjects related to interaction between
adsorbed atoms and surfaces. They are the Atomic Force Microscopy and atomic-

scale friction at graphite surface, and the structure and dynamics of Pd-H systems.

In Chapter 2, I will review the computational tools I used. I will start with the
Density Functional Theory, in particular its implementation in the Local Density
Approximation (LDA) technique. I will briefly discuss the computational details,
such as the basis and use of ab initio pseudopotentials. Then, I will derive the Many-
Body Alloy Hamiltonian for metal and alloy systems. At the end, I will present the
formalism for Nosé and Rahman—Parrinello Molecular Dynamics simulations, which
allows for the description of the dynamics of a canonical ensemble exposed to an

external anisotropic stress field.

In Chapter 3, I will develop the theory of Atomic Force Microscopy. I will first
present the ab initio results for the Pd-graphite interaction. Next, I will present a
parametrization of the Pd-graphite interaction potential, which is inspired by Embed-
ded Atom Method and which is not restricted to high symmetry sites. Combined with
continuum elasticity theory, these results will be used to determine whether atomic
resolution can be achieved nondestructively in AF M experiments on the graphite sur-
face. Then, I calculate the atomic scale friction of a sharp Pd tip on graphite and
determine the friction coefficient. I propose two idealized friction machines to explain

the possible microscopic mechanisms of energy dissipation in the friction process.

In Chapter 4, I calculate the equilibrium structure and dynamics of Pd-H systems.

This chapter is divided into two parts. In the first part, I construct and test the many-

12

body alloy Hamiltonian for the Pd-H system. I calculate the equilibrium structure
and surface phonon spectra of clean and hydrogen covered Pd surfaces, in order to
determine the effect of hydrogen on the equilibrium structure and dynamics of Pd.
In the second part, I use molecular dynamics to study the dynamical properties of
Pd-H systems at finite-temperature. I will study the melting transition of pure Pd in
the bulk bulk and the mechanical stability of bulk Pd at different temperatures and
hydrogen concentrations. I will show that the mechanisms of hydrogen embrittlement

can be understood by carefully analyzing the MD simulation results.

Chapter 2

Theory

2.1 Ab initio Density Functional Formalism

The basic theory governing the electronic and structural properties of solids is quan-
tum mechanics. In systems with many particles (such as ~ 1023 in solids), the exact
solution of quantum mechanical equation is essentially impossible to obtain, due to
the complex many-body interactions which couple many degrees of freedom. Even
though the ground state energy can be expressed in terms of one- and two-particle
Green’s functions, the computation of these quantities involves a set of differential
equations which couple all n-particle Green’s functions. The basic difficulty to de-
scribe even ground-state properties of solids is significantly reduced in the Density

Functional Theory.

2.1.1 Density Functional Theory

In this Thesis, I will derive the Density Functional Theory (DFT) formalism very
briefly, and refer the reader to other excellent reviews for more details [Lun 83,
Cal 84, Jon 89, Mah 90, Dre 90]. The DFT is based on two theorems proposed by
Hohenberg and Kohn in 1964 [Hoh 64], which state that:

13

14

l. The electron density n(1"') in the ground state is a functional of the potential

no

2. The potential V0") is a unique functional (to within a constant) of the electron

density n(r').

This is equivalent to saying that the exact ground state properties of a system can
be calculated using a variational approach involving only the electron density, rather

than an antisymmetric wave function.

I follow the derivation by Levy [Lev 79]. The Hamiltonian describing the motion

of N electrons in an external potential me is
H=T+%+Km on

where T and Vcc are the kinetic and electron-electron interaction operators, respec-
tively. For all densities n(r') which can be obtained from an antisymmetric wave

function 112(r'i, r3, . . . , rTv), Levy defined the functional
Fin] = minWIT + Veelib) = ( 3.;an + Veell‘pihin), (2-2)
where the minimum is taken over all 1b that give the density n, and tbmgn minimize

F[n].

Let us denote Egg, $35, and nas(r‘) to be ground state energy, wave function, and

electron density of the system under external field Vm. Then,

Bin] 2 [dv’Vmb’inb’HFin] (2.3)

( 3“an + We ‘1" I/extlfpr’hin)

> 5'05, (2-4)

.

15

which follows from the variation principle applied to the ground state. This relation

applies in particular to the ground state wave function swhich gives

WIGSIT + Vcc + thl'l’cs) S (1033.?ng + Vise + Kali/1333:)- (25)

"GS

The charge density corresponding to $35 and lbw-n are the same, so the interaction

with the external potential can be subtracted on both sides of the inequality,

WGSIT + Welt/’05)- W) 33373” + Kali/’i'lflrfl- ‘ (2-6)

Combined with Eq. (2.4), we find

(thoslT + Kelwcsl= < 333.? IT + Welt/’33:) = Fines} (27)

So, the ground state energy

Ens = fdmzdflnosbf) + F[nos] (2.8)

is a functional of ground state charge density. This equation, combined with Eq. (2.4)
completes the proof of the basic theorems. These theorems provide a general method
for calculating ground state properties. However, theses theorems do not suggest

any particular form of the functional F [n] To solve this problem, Kohn and Sham

introduced the Local Density Approximation (LDA) in 1965 [Koh 65].

2.1.2 Local Density Approximation

In their famous paper [Koh 65], Kohn and Sham decomposed the energy functional

as

E[n] = To[n] + F + E“, (2.9)

16

where

(7'1") (7‘2)

[(7‘1 - r2].

2W. 1")n( 1")+//dr'idr'§n (2.10)

Here, To is the kinetic energy of a system of non-interacting quasi-electrons with the
density n(F). F is the Coulomb energy of the system of quasi-electrons. E“ is the
exchange-correlation energy, which gives the contributions to the total energy except

for the kinetic and Coulomb terms.

Variational theory, applied on this energy functional, gives

T__o__[n]6

Tn +I/ezt‘l'VH'f'ch—I-‘zo (211)

Here, V” is Hartree potential derived from the second term of F, VIC = 6Exc[n]/6n,
and u is the Lagrange multiplier which constrains the total number of particles to be
constant. On the other hand, the variational equation for a system of non-interacting

quasielectrons reads

6To[n] + V

6n - p = 0. (2.12)

Inspection of Eqs. (2.11) and (2.12) shows that they are the same if V: chg+VH+Vw
The solution of Eq. (2.12) can be obtained using a Schr6dinger equation. For a system

of quasi-electrons, this leads to a set of Kohn-Sham equations
R2
l-gng-V’ + v... + vim + molt-(o = 6.1M") (2.13)
and
"(a = z [WU-"llz- (2.14)

Note that e,- and 1b, are not necessarily the eigenvalues and eigenfunctions of the

system; only the electron density n(1'~') is the correct physical quantity.

17

Kohn-Sham equations involve a universal functional, the exchange-correlation en-
ergy E,c[n], and its functional derivative V1.60"), which are only known exactly for a
homogeneous electron gas. As a manageable simplification which has been introduced
in the Local Density Approximation (LDA), the functional Erc[n] is replaced by a

local function Exc(n). Hence,

E“ = fdf'negc(n(f")), (2.15)

where age is the exchange-correlation energy density for the homogeneous electron gas.
Many parametrizations of file have been proposed, such as that of Wigner[Wig 38],
Kohn and Sham[Koh 65], Hedin and Lundqvist[Hed 71], Ceperly and Alder[Cep 80]
etc.. In our calculation, we usually use Hedin-Lundqvist parametrization, which have

been proven to give good results.

The LDA is a accurate in systems with a slowly varying charge density. It is also
accurate in systems with a large charge density, where the kinetic and Hartree terms
dominate E“, and in systems with small charge densities, where gradient of n(i~')
are small. This formalism has been successfully applied to many systems, including

atoms, molecules, clusters, surfaces and bulk of solids.

2.1.3 Computational details

In the real calculations based on the DFT, the initial guess of the charge density is
obtained using a superposition of atomic charge densities. Using this charge density as
input, V“ can be obtained using the LDA. This provides all of the information needed
to set up the Kohn-Sham equations, which are usually solved as matrix equations in
a given basis. The solution of the Kohn-Sham equations gives e,- and 11),. Then the

Fermi energy can be determined and the charge density obtained as output. In the

18

next iteration, this charge density is used as input, and the procedure is repeated
until self-consistency is achieved. In the self-consistent field method, the input and

output potential is required to be the same.

One practical consideration in a realistic calculation is the selection of a basis. The
plane wave basis is the simplest and easiest to implement, since all basis functions
are orthogonal to each other. Plane waves give a complete set of basis functions;
a finite basis is typically limited by the cutoff energy EC. This basis is ideal for
nearly free electron systems like simple metals, where the effective potential does
not vary too much. The number of basis functions needed is proportional to Bf”,
where d is the dimensionality of the system. The computational load increases as
Eff/2, since solution of Kohn-Sham equations involves matrix diagonalization. This
computational effort makes a plane wave basis impractical for materials with localized

electrons, i.e. transition metals and semiconductors.

In my Linear Combination of Atomic Orbitals (LCAO) [Laf 71, Cal 72, Che 84]
calculation, I use a local Gaussian basis in order to improve efficiency [Har 82]. The

form of the basis functions is

ezp(—a,-r2) for 3 states

f.- = (2:,y, z)ea:p(—a,-r2) for p states

(:ry, yz, 22:, a:2 — y2, r2 — 322)e.rp(—a,-r2) for d states
This basis provides an efficient representation of the atomic wave functions in
systems with localized electrons. More than one radial Gaussians are typically used
for each orbital to allow for variational freedom. Completeness of the basis is tested
by increasing the number and location of Gaussian functions and monitoring the
convergence of the calculated total energy. In my calculations, usually three to four
Gaussian decays for each atomic orbital (s, p,, py, ...) are found to give a sufficiently

complete basis set. The decay constants a for a given atomic type are chosen to

 

19

minimize the total energy of the bulk system. Since the LCAO basis contains non-

orthogonal basis functions, the solution of the Kohn-Sham equations is more involved.

Beside using Gaussian functions for basis functions, it is very useful to fit the
screened atomic potentials and charge density to a set of Gaussians as well. This
procedure allows for an efficient evaluation of crystal LDA-Bloch functions up to very

high Fourier components corresponding to a large value of the energy cutoff EC.

Most electronic and structural properties of solids are determined by the behavior
of valence electrons, while the core electrons are essentially inert. Valence electrons
feels a very weak potential outside the core due to the screening by the core electrons.
Replacing the true ionic potential by a pseudopotential, which models the nuclear and
the nonlocal core potential, simplifies the calculation considerably. I use the ab initio
pseudopotential form proposed by Hamann, Schliiter, and Chiang [Ham 79, Bac 82],

which has following desirable properties:

1. Pseudopotentials are continuous non-divergency functions with a continuous

first derivative everywhere.

2. Real and pseudo eigenvalues for valence states agree for a chosen “prototype”

atomic configuration.
3. Pseudo wavefunctions are nodeless.
4. Real and pseudo wavefunctions agree beyond a chosen “core radius” re.

5. The integral from 0 to r of the real and pseudo charge densities agree for r > rc

for each valence state (norm conservation).

6. The logarithmic derivatives of the real and pseudo wave functions and their first

energy derivatives agree for r > rc.

20

Properties 5 and 6 are important to ensure optimum t‘ransferability of the pseudopo-
tential between different chemical environments. Norm conservation enables the use
of the pseudo charge density instead of the real charge density in UP T calculations.
Property 6 ensures the correct scattering properties of the potential. No adjustable
parameters occur in the construction of the Hamann-Schliiter-Chiang pseudopoten-

tials.

The LCAO Gaussian basis and ab initio pseudopotentials are very powerful which
make precise LDA calculations for transition metals or semiconductors feasible. Nev-
ertheless, the computational load increases rapidly with increasing complexity of the
system (reduced symmetry, increased size of the unit cell). For very large systems,
parametrized calculation schemes are more appropriate. A very efficient scheme to
calculate the total energy of a system is the Many-Body Alloy Hamiltonian to be

discussed below.

2.2 Many-Body Alloy Hamiltonian

The Many-Body Alloy (MBA) Hamiltonian is an extension of a total energy scheme,
which has been successfully used previously to study the electronic and structural
properties of small clusters, surfaces of metals and dilute metal alloys [Tom 83,
Tom 86a, Tom 85b, Tom 85a, Spa 84]. Its derivation has been published recently
by W. Zhong, Y.S. Li, and D. Tomanek [Zho 91a].

As discussed earlier [Tom 83, Tom 86a, Tom 85b, Tom 85a], the total cohesive
energy of a solid can be decomposed into individual atomic binding energies Ecoh(i),

8.8

E...<tot) = 2: use). (2.16)

21

The binding energy of atom i consists of an attractive part due to the hybridization

of orbitals, E33 (i), and a term ER(i) describing repulsive interactions. Hence,
' E,,,.(i) = 1335(2) + 133(2). (2.17)

Different simplified parametrization forms have been proposed [Bre 89] for the many—
body energy EBS(i). The embedded-atom method (EAM) [Daw 84, Foi 85, Foi 87,
Daw 89, Foi 89] takes EBS(i) as a unique function of the total charge density of the
unperturbed host at the site i. Since this parametrization might cause problems in
the case of alloys with nonzero charge transfer, I base my expression for E85 (i ) on a
tight-binding Hamiltonian. In a one-electron picture, the binding energy of atom i is

given by an integral over the local density of states at i, N,( E), as
Er
535(2) = — j (E - E0)N,-(E)dE. (2.18)

In the second moment approximation, EBS(i) is proportional to the effective band-
width, which in turn is proportional to the square root of the second moment M2(i)
of the local density of states. Then,
1/2 . 1/2
E35(i) or M5(i)1/2 = {2%.} oc (Zia-2W} . (2.19)
1'?“ is“

In the last part of this equation, I have related M2 to the hopping integral tij between
neighboring sites i and j and assumed an exponential distance dependence of the

effective (screened) hopping integrals, as t(r) o< e’".

The repulsive part is parametrized by a pairwise Born-Mayer potential with an

exponential distance dependence, as

153(2) or 2*: (W. (2.20)

22
Then,
13.0,.(2) .—. E35(i)+ER(i)

n.1, 1/2
-(253 expl-2q(;; - 1M} +65 2 expl-p(

j?“ jgéi

7‘

r“ —1)]. (2.21)

Here, rgj is the distance between atoms i and j. Parameters p and q describe the dis-
tance dependence of the hopping integrals and the Born-Mayer interactions, respec-
tively, and are related to the bulk elastic properties. In the case of a single-component
bulk crystal, Eq. (2.21) can be used to reproduce the equilibrium properties, such as
the equilibrium nearest neighbor distance r0 and the bulk cohesive energy Ecoh(bulk).

In this case, assuming isotropic hopping integrals, £0 and E32 in Eq. (2.21) are given

 

by
_ Eco),(bulk)
£0 _ (1 -<I/10)(Zbun:)‘/2 ’ (2'22)
65 = 5° 1, (2.23)

 

2.11;, p
where Zhu“, is the bulk coordination number.

For systems with more than one component (such as alloys and compounds),

Eq. (2.21) can be generalized to

1" 1/2
Ecoh(iy 0) = -{Zj¢i €3,afi eXpl—2qafi(-u££ — 1”}

70.05

+ 21%|. 661,03 8Xp[—pag(m _ 1)]7 (2.24)

'0,afi
where a, 6 represent the types of atoms i and j, respectively.

As I discussed above, the Hamiltonian underlying the energy expression in
Eq. (2.24) describes the essential physics governing cohesion in many solids. The
large flexibility and the microscopic basis for the description of many-body attractive

interactions in alloys makes the MBA Hamiltonian superior to embedded-atom like

23

schemes [Daw 84, Bre 89] for two main reasons. First, unlike the EAM, Eq. (2.24)
does distinguish the changed binding of atom i that is surrounded either by Z, atoms
of type a or Z3 atoms of type 6. This holds even in the case that the charge density
at site i due to the surrounding atoms is the same. Second, in contrast to the EAM,
this approach does not assume local charge neutrality in alloys. For this reason, I
feel confident to apply this energy expression as an “intelligent interpolation scheme”
to determine the energy of structures with low symmetry, once the corresponding
ab initio data for high symmetry structures are available. Specifically, I will use
this energy expression to study the equilibrium structure and dynamics of selected
systems. The calculation of dynamical properties will focus on phonon spectra, to
be addressed by lattice dynamics techniques, and large scale deformations during the

process of melting and fracture to be addressed in molecular dynamics simulations.

2.3 Lattice dynamics

Lattice dynamics is used to determine phonon dispersion relations of a crystaline
solid. The phonon spectra can be compared directly to experimental results for the
equilibrium state at very low temperatures. The phonon frequencies of a crystal can
be obtained using the dynamical matrix [Mar 63]. This matrix D(li), corresponding

to the wave vector If, is given by

D..,,,,..,(i£) = (M.M.,)-1/2 Z oafi,,.,,-.,exp[—ii’.(it — 1%)]. (2.25)
a—a

Here, M“ is the mass of the n-th atom and My is the mass of the u-th atom. $0.53,,“

is the force-constant matrix, which can be expressed as

 

Q . . _ B’Ew).(tot)
aflnngu — aua,5,‘6Ufi,jy a

(2.26)

24

where Ecoh(tot) is the total energy of the system. um,“ is the a-th Cartesian compo-
nent of the displacement of n-th atom in i-th unit cell, and Ugh“, is the fl-th Cartesian

component of the displacement of the u-th atom in the j-th unit cell.

Finally, the phonon frequencies w(li) are the eigenvalues of DUE), which are cal-

culated by solving

det|w2(li)l — D(i£)| = 0. (2.27)
2.4 Molecular dynamics

Unlike lattice dynamics, molecular dynamics (MD) describes the evolution of the sys-
tem with time, based on a direct numerical solution of the equations of motion for
individual atoms. Atomic motion is not restricted to the close vicinity of the equilib-
rium site. Such a simulation provides microscopic information about the dynamics of
the system, including structural changes during phase transitions (such as melting or

stability breakdown under tensile stress).

The dynamics of an isolated N -particle system in three-dimensional space is gov-

erned by the Lagrangian

N
1 .

L = 2 5mm? - V({Qi}) - (228)

i=1
Here, q is the position vector of atom i in the system and V({q,}) is the total poten-
tial energy, given by the MBA Hamiltonian. The dynamical evolution of the system
is described by Euler-Lagrange equations derived from the above Lagrangian. This
procedure yields statistics for a micro-canonical ensemble. Unfortunately, the simu-

lation of realistic micro-canonical systems requires a very large number of particles

which imposes unrealistic requirements on the computation.

25

A more natural choice is the canonical ensemble, which considers the temperature
(rather than the total energy) to be constant, and which allows for a free energy
exchange with the external heat bath. I use an algorithm due to Nosé [Nos 84] to
describe a canonical ensemble. In the N osé scenario, the fixed temperature canonical
ensemble can be simulated by considering a single additional variable 3. In this
extended phase space, the net effect of the heat bath is assumed to be a scaling of
the velocities as vi = sqj. Nosé interpreted the scaled velocity V, as the true velocity
resulting from the heat exchange with the external heat bath. The canonical ensemble

can be described in an augmented Lagrangian of the form

N 1 . 1 , '
L = E 5"2,92qu — V({qj}) + 503’ — (N, +1)T,,1ns . (2.29)

i=1
Here, N f = 3N is the number of degrees of freedom of the system. Q is the mass
associated with the new variable 3 which describes the coupling to the external heat
bath. T5,, is the temperature of the heat bath and consequently the equilibrium
temperature of the system. The optimum choice of Q provides an efficient damping
of the equilibrium state. Overdamping is avoided by choosing Q in such a way
that the oscillation period of s is much smaller than the damping time constant.
The Lagrangian in Eq. (2.29) contains a logarithm of s in the potential energy part
which makes this extended system equivalent to the canonical ensemble describing the
original N -particle system. This new Lagrangian is justified by the correct statistics
it yields for the system. In other words, when integrating the partition function over
the new dimension in phase space, I regain the partition function for the canonical
ensemble. Consequently, the microcanonical ensemble of the augmented Lagrangian,
described by Eq. (2.29), generates precisely the canonical ensemble of the original

N -particle system.

In order to describe the effect of mechanical coupling between the ensemble and

26

the outside world, Andersen applied molecular dynamics to systems under constant
hydrostatic pressure [And 80]. In this technique, the variable space is extended by
an additional quantity which scales all the atomic coordinates in the unit cell, thus
allowing the volume to change uniformly. This method can be combined with Nosé’s
treatment of the canonical (TN) ensemble to describe the behavior of that system at
constant tensile stress (or tension) t (TtN ensemble). This method has been further
generalized by Parrinello and Rahman [Par 80, Par 81, Par 82, Ray 84, Ray 88] to

allow for anisotropic stress and corresponding shape changes of the unit cell.

The mechanical stability of Pd and PdH systems which will be addressed in Chap—
ter 4 is intimately related to the dynamics of these systems under uniaxial tensile
stress. I will address this process using the method of Parrinello and Rahman [Par 80].
In the molecular dynamics calculations, I will consider a unit cell which is spanned
by the three Bravais lattice vectors a, b and c. These vectors define a shape matrix
11 = (a, b, c) which is used to describe the response of the system to applied stress.
This matrix maps the true atomic coordinates onto reduced coordinates q; which lie
in a cube of unit length, so that the true coordinates are obtained by the scaling

transformation hqj. The Lagrangian describing the TtN ensemble is then given by

N

L = ng.s’(qiqui)-V({q1})+gas—(NH911.1113
i=1
+ §W.Tr(hfh)—%Tr(u). (2.30)

The first term describes the kinetic energy of the system, and the matrix Q is given
by Q = hTh. W5 is the mass of the “piston” which exerts the external tension on the
system. The equilibrium volume of the MD unit cell at zero tension is denoted by

V0 and the equilibrium shape by 110- In the case of nonzero tension, I define a strain

27

matrix g by

[(ho-1)T§ho'l - I] . (2-31)

NIH

g:

The quantity 1 in Eq. (2.30) is the thermodynamic tension tensor which gives the
work dW due to an infinitesimal distortion of the system as dW = %Tr(t_e). t is

related to the stress tensor g by [Par 80]

a = %h(ho“).t.(ho“)ThT/V. (2.32)

where h is the average value of the shape matrix and V = det(h) is the average

volume.

Once the Lagrangian L is established, the equations of motion are given by

_¢_1_ 6L _ aL
dt 345a aqia
i i _ PL _. 0
dt 81:0,, 340:9 -

2.2912. -0
dtas as_'

0

 

 

(2.33)

Here, I have used Greek indices to denote the vector and matrix components of Qj

and h, respectively.

The dynamics of the system is obtained by integrating the equations of motion in
real time using a predictor-corrector method of fifth order [Gea ]. The typical time
step in the MD calculation of bulk Pd is 2 x 10‘153. For hydrogen loaded Pd, the
time step is reduced to 5 x 10‘163 due to the high vibration frequency and large-scale
diffusion of the light hydrogen atoms. The generalized energy of the system is a

constant of motion. I control the numerical fluctuations of this quantity (due to the

28

finite time steps) to lie within a small error margin of 10", which is indicative of
a precise integration. In order to avoid overdamping, I tune the mass Q associated
with the temperature and the mass Wt associated with the tensile stress in such
a way that the fluctuations of temperature and volume occur over large periods of
2150-100 time steps. In a typical MD calculation, I consider a periodic arrangement
of unit cells (with originally cubic shape) which contain 500 Pd atoms and a fixed
number of hydrogen atoms. When studying the properties of the system at a specific
temperature or pressure, I first let the system equilibrate over a period of z30,000
time steps and only then start collecting data for statistics. I find this time period to

be sufficient for an excellent statistics with negligible error bars.

In the derivation of the algorithm, Nosé has assumed the system to be ergodic
[Nos 84]. Doubts have been raised about the validity of this assumption [Kus 90], es-
pecially in the case of small systems. There are two reasons why this criticism should
not affect the results presented in Chapter 4.3. First, I consider a very large system
with more than 1500 degrees of freedom. Second, the many-body alloy potential con-
tains anharmonic terms which gain importance at interatomic distances substantially
different from r0. Especially at elevated temperatures, this fact contributes to a fast

onset of chaotic behavior.

Chapter 3

Theory for the Atomic Force
Microscopy

As discussed in Chapter 1, the Atomic Force Microscope (AFM) is a recently devel-
oped instrument which allows a direct imaging of topological structure at insulator,
semiconductor, and metal surface alike [Bin 86, Bin 87]. The challenging question is
whether the AF M can achieve atomic resolution of a given surface, and what are the
theoretical limitations of this capacity. The study of the interaction between an AFM
tip and a surface are also of relevance for the fundamental understanding of friction

and the related phenomena.

These problem areas which are the subject of this Chapter, have been addressed

in the following seven publications:

1. W. Zhong and D. Tomanek, F irst-Principles Theory of Atomic-Scale Friction,
Phys. Rev. Lett. 64, 3054 (1990).

2. W. Zhong, G. Overney, and D. Tomanek, Theory of Atomic Force Microscopy on
Elastic Surfaces, The Structure of Surfaces III, edited by S.Y. Tong, M.A. Van
Hove, X. Xide and K. Takayanagi, Springer-Verlag, Berlin (1991), p243.

3. W. Zhong, G. Overney, and D. Tomanek, Limits of Resolution in Atomic Force
29

30

Microscopy Images of Graphite, Europhys. Lett.‘ 15, 49 (1991).

4. G. Overney, W. Zhong, and D. Tomanek, Theory of Elastic Tip-Surface Inter-

actions in Atomic Force Microscopy, J. Vac. Sci. Technol. B9, 479 (1991).

5. D. Tomanek, W. Zhong, and H. Thomas, Calculation of an Atomically Mod-
ulated Friction Force in Atomic Force Microscopy, Europhys. Lett. 15, 887
(1991).

6. D. Tomanek and W. Zhong, Palladium-Graphite Interaction Potentials Based

on First-Principles Calculations, Phys. Rev. B43, 12623 (1991).

7. G. Overney, D. Tomanek, W. Zhong, Z. Sun, H. Miyazaki, S. D. Mahanti, and
H.-J. Gfintherodt, Theory for the Atomic Force Microscopy of Layered Elastic
Surfaces, J. Phys. : Cond. Mat. 4, 4233 (1992).

The description of tip-surface interactions in the AFM is determined quantitatively
for Pd/ graphite as a prototype system. Pd is a non-magnetic transition metal with
near-noble electronic configuration. The valence 4d electrons are less localized than in
3d elements, yet unlike in the 5d elements, no relativistic effects need to be considered.
This element is also chosen for its wide spread industry application. Graphite is a very
common and important material in the industry and research alike. The inert nature
of graphite layers facilitates cleavage and preparation of perfect surfaces. Graphite
is also one of the most commonly used systems in AFM and STM studies of surface

morphology [Mey 88].

31

3.1 Calculation of the Pd-graphite interaction

Accurate determination of the interaction energy between Pd and graphite is crucial
for predictive calculations of AFM images of graphite, obtained with a Pd AFM tip.
I start with an ab initio calculation of the interaction between a perfect layer of Pd
atoms and the graphite substrate. In this calculation, I determine the adsorption
energy of Pd atoms in the on-top (T) and the sixfold hollow (H) site as a function
of the Pd-graphite separation. Total energies are calculated from first principles
using the Local Density Approximation (LDA) [Koh 65]. I am using the ab initio
Pseudopotential Local Orbital Method which has been described in Chapter 2 and
successfully applied to short wavelength distortions of graphite [Tom 89]. In this
calculation, ionic potentials are replaced by norm-conserving ionic pseudopotentials
of Hamann—Schlfiter-Chiang type [Ham 79] and the Hedin-Lundqvist [Hed 71] form
of the exchange-correlation potential is used. The surface of hexagonal graphite is
represented by a 4—layer slab and the adsorbate by a monolayer of Pd atoms in registry
with the substrate (1 Pd atom per surface Wigner-Seitz cell of graphite). The basis
consists of s, p and d orbitals with four radial Gaussian decays each on Pd sites and
of s and p orbitals with three radial decays each on carbon sites, i.e. 40 independent
basis functions for Pd and 12 basis functions for C. In order to obtain accurate energy
differences between the “T” and “H” geometry, I introduce floating orbitals on sites
not occupied by atoms and use the same extended basis [Tom 86c] for both “T” and
“H” calculations. To insure high accuracy and complete convergence of total energies,
I use an energy cutoff of 49 Ry in the Fourier expansion of the charge density and
sample the 2-dimensional Brillouin zone with a mesh of 49 I'd-points, using a special-

point scheme [Cha 73].

In Fig. 3.1, I show the adsorption energy Ead of a Pd atom (representing the Pd

32

 

 

 

 

Ead(ev)
HOHNUOFCJ‘OJ'QCDCDO
l

 

 

 

 

 

 

1.0 , 1.5 2.0 2.5 3.0 3.5

Figure 3.1: Pd adsorption energy E“ as a function of the adsorption height 2 above
the surface of hexagonal graphite. The solid and dashed lines correspond to the
sixfold hollow (H) and the on-top (T) sites, respectively. An enlarged section of the
graph near equilibrium adsorption is shown in an inset. A second inset shows the
adsorption geometry in top view; a possible trajectory of the Pd layer along :1: is
shown by arrows. (From Ref. [Zho 90]. ©American Physical Society)

33

monolayer), defined by
E...) = Etom(Pd/graphite) - Etom(Pd) - E,o,d(graphite),

for the on-top and hollow site registry. The calculation yields a very weak ad-
sorption bond strength 30.1 eV at an equilibrium adsorption height 2 z 3 A. At
this distance from the substrate, the corrugation of the graphite charge density is
negligibly small due to Smoluchowski smoothing and the position-dependence of the
adsorption bond strength is << 0.1 eV. This effect is also responsible for the large
value of the surface diffusion constant and apparent small sensitivity of surface fric-
tion to adsorbed films as discussed in Ref. [Ski 71]. These calculations show that
at bond lengths 2 S 2 A, the hollow site is favored with respect to the on-top site.
At 2 z 2 A, the adsorption energies are nearly the same, and at larger distances,
it is the on-top site which is slightly favored by < 0.05 eV. This is consistent with
the dominant interaction changing from closed-shell repulsion (which strongly favors
the hollow site at very small adsorption bond lengths) to a weak chemisorption bond
(which is stabilized by the hybridization with p, orbitals in the on-top site). The
magnitude and site dependence of van der Waals interactions between AF M tip and
the surface is negligible in the weakly repulsive region of the Pd-graphite potential

considered here.

The DFT calculations provide us with very accurate results, but only for small
systems with high symmetry. When the Pd atoms are located at low symmetry
sites, the calculation is computationally much more intensive. In calculations of
AFM images and atomic-scale friction, more efficient total energy schemes have to be
used. Model potentials for the interaction with graphite, the material of interest in
this study, are typically based on interactions with a “generic” carbon solid [Gou 89,

Abr 89]. These model potentials clearly cannot distinguish the surface reactivity of

34

diamond with an sp3 configuration from that of graphite with sp2 bonding.

Here, I present a simple expression for the metal-graphite interaction potential
which is based on a first-principles total energy calculation of Pd on graphite [Zho 90].
While derived specifically for the Pd/ graphite system, the form of the interaction
potential is more general and is a good prototype for the interaction of any metal
adsorbate on graphite. The expression for the metal-graphite interaction, inspired by
the Local Density Approximation (LDA) [Koh 65], provides a basically correct picture
of many-body interactions in the adsorption system. The simple parametrization is
a major advantage which will allow this potential to be used in computationally

intensive molecular dynamics simulations.

It is convenient to note that the LDA results for Edd can be well approximated by
a local function which depends only on the total charge density of the graphite host

at the Pd adsorption site,

EMU") = Ead(P(O)- (3-1)

This form of the interaction potential is inspired by the the Density Functional for-
malism [Koh 65] and the Embedded Atom Method [Daw 84, Foi 85, Foi 87, Daw 89,

Foi 89]. Then Ead(p(i")) can be conveniently parametrized as

EMPW) = 51 (p/po)°" — 62 (P/P0)a’- (3-2)

In the case of Pd on graphite, £1 = 343.076 eV, 52 = 2.1554 eV, a1 = 1.245, 02
= 0.41806, and p0 = 1.0 e/A3. The dependence of E“) on p, obtained using the

parametrized form in Eq. (3.2), is shown in Fig. 3.2.

From an independent LDA calculation of graphite surfaces, I find that the total

35

 

0.5
40 " 0.7

35 " 0°” , 1

I
J

U

 

End (9")

 

 

 

 

 

 

 

 

L L L L l L

1 2 3 4 5 6 7 6 910

-1n(p/el.§'°)

 

Figure 3.2: Relation between the Pd adsorption energy 505(5) and the total charge
density of graphite p(r"') at the adsorption site F, given by Eq. (3.2). An en-
larged section of the graph near equilibrium adsorption is given in the inset. (From

Ref. [Tom 91b]. ©American Physical Society)

36

charge density can be well approximated by a superposition of atomic charge densities,

pm = 25,4;— 1‘2’.) (33)

This parametrization is especially convenient in the case of deformed surfaces
where an LDA calculation is difficult due to reduced symmetry. On flat surfaces, the
maximum difference between the LDA charge density and the superposition of atomic

charge densities is only few percent.

Finally, it is useful to parametrize the charge density of carbon atoms. I find that

pat, as obtained from an atomic LDA calculation, can be conveniently expressed as
pai(r) = pce’B', (3-4)

where pa = 6.0735 e/A3 and fl = 3.459 A”. As shown in Fig. 3.3, this expression is
a very good approximation to the LDA results especially in the physically interesting

range 1.4 A < r < 3.0 A.

One of the primary uses of the above potential is to describe the interaction
between a metal AF M tip and graphite. In case of a large tip, one can expect a
contribution to the interaction potential from Van der Waals forces, which are not
described correctly by LDA especially at large tip-substrate distances. As discussed in
Chapter 3.2 and in Ref. [Zho 91b], these forces are not very important since they are
very small (typically < 10'10 N) at tip-substrate separations 2,3 A. At smaller separa-
tions, they are dominated by the substrate-tip repulsion which is described correctly
within LDA. These dispersive forces are also not very important for the interpretation
of experimental results, since they do not show atomic resolution [Hor 92] and are
easily compensated in the experiment by adjusting the force on the cantilever which

supports the tip.

37

 

—ln(p/el.fi’°)

 

 

 

 

2 L L L l L L
-

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5

r (A)

Figure 3.3: Radial plot of the charge density of a carbon atom p550"), based on an LDA
calculation (dashed line). The solid line shows the parametrized form of the charge
density, given by Eq. (3.4). ( From Ref. [Tom 91b]. ©American Physical Society)

38

A parametrized form of the potential, given in Eqs. (3.2), (3.3) and (3.4), allows for
a very fast and efficient evaluation on a computer and is ideally suited for molecular
dynamics calculations. Adsorption energies of Pd on a flat graphite surface, based on
this parametrized form, are compared to LDA results in Fig. 3.4. I find a very good
agreement between LDA results and parametrized potential over the whole interaction
range. The largest discrepancy occurs at z = 2.0A, where the switching between the
preferential on-top and hollow sites is not reproduced by the parametrized potential,
corresponding to an error of S 0.1 eV. In the following, I will use the ab initio results

whenever possible, and the parametrized results in other cases.

Summarizing the results of this section, I calculated the Pd/ graphite adsorption
energy as a function of the adsorption height, for Pd at the on-top and hollow site.
The interaction is characterized by a hard-core repulsion at small distances and very
weak chemical bond at larger distance (2 > 2.5A). This interaction can be well

parametrized in terms of the charge density of the graphite host.

3.2 Limits of resolution in AFM images of gra-
phite

As pointed out in Chapter 1, the possibility of atomic resolution in the AFM is
limited in two ways. Too small forces lead to insufficient corrugations, and too large
forces may destroy the surface. So far, these limitations of the atomic resolution have
not been addressed in the literature. Present theoretical information is limited to
calculations of the interaction between an infinite “periodic” carbon or aluminum tip
and a rigid surface [Cir 90, Bat 88], and the interaction between a single AFM tip and
an elastic surface represented by a semi-infinite continuum [Tom 89] or by a model

system of finite thickness [Abr 89, Lan 90].

39

 

End (9")

 

 

 

 

 

 

 

 

Figure 3.4: Adsorption energy E“ of Pd on graphite as a function of the height 2 of
Pd atoms above the graphite surface. F irst-principles results are given by 0 and e for
the hollow and the on-top site, respectively. These data are compared to the present
results, obtained using Eqs. (3.2), (3.3) and (3.4) and given by the solid and the
dashed line for the hollow and the on-top site, respectively. Upper inset: Schematic
top view of the adsorption geometry. Lower inset: An enlarged section of the graph
near equilibrium adsorption. (From Ref. [Tom 91b]. ©American Physical Society)

 

40

In this Section, I discuss for the first time the theoretical limits of atomic resolution
in AFM. I will use the ab initio results of the interaction between a “sharp” monatomic
Pd tip and the surface of graphite discussed in Section 3.1, to predict corrugations
for a varying load Fm applied on the AFM tip. These calculations also predict
tip-induced elastic substrate deformations which limit the range of applicable loads

Fert-

Using the adsorption energies End given by an ab initio calculation, the force on

a “sharp” monatomic AFM tip is given by

_ aEad(z)
fez-t — '— 62 ~ ' (3'5)

 

Here, microscopic forces per atom are denoted by f and macroscopic forces applied
to large objects by F. I also investigated the effect of the long-range Van der Waals
forces on the tip-substrate interactions which are not described correctly by the LDA

especially at large tip-substrate distances.

The Van der Waals force between an extended conical tip and a flat surface is
estimated using the expression Fv4w(z) = A3 x tanza/(6z), where a is half the
opening angle of the cone [And]. In this expression, A” is the Hamaker constant
and z is the distance between the cone tip and the surface. In this calculation, I
consider a = 30° and A3 = 3 x 10‘19 J, which is a typical value for metallic systems.
For tip-substrate distances 2 > 3 A, the Van der Waals forces are very small, typically
Fww < 10-10 N. At smaller distances, these forces can be neglected when compared
to the closed-shell and internuclear repulsion which are described correctly within
LDA. Since each of these regions is dominated only by one type of interaction, I
determine the total tip-substrate force Fm as a superposition of the force described

by LDA and the Van der Waals force. In Fig. 3.5(a), I present results for the total

41

force F = fez, between a conical AFM tip with 1 Pd apex atom and a graphite

substrate, together with a schematic top view of the geometry.

The LDA calculations yield nearly the same weak Pd-graphite interaction poten-
tial (Edd < 0.1 eV) in the “H” and “T” sites beyond the equilibrium Pd-graphite bond
length 2,, z 3 A and consequently the same weak interaction force. This is plau-
sible since End is closely related to the total charge density )0, shown in Fig. 3.5(b),
which has only a very small site-dependence in the attractive region of the potential
at z > 2,, due to Smoluchowski smoothing [Smo 41]. Consequently, the corrugation
A2 in this region is too small to be detected by the AF M. In the weakly repulsive
region of the potential, for 2 A < z < 3 A, the on-top site is slightly favored with
respect to the hollow site due to a weak chemisorption bond with substrate C2125
orbitals. For atom bond length 2 < 2 A, the strongly repulsive Pd-graphite interac-
tion is mainly determined by the closed-shell repulsion which energetically favors the
hollow site. Hence, in the repulsive region of the potential, the Pd tip comes closest
to the substrate near the “T” site for small loads. For large loads, the tip is closest

to the surface near the “H” site.

Fig. 3.6(a) shows the expected AFM corrugation Az during an my scan of the
graphite surface, for feet = 10‘8 N. The equilibrium bond length 2(fm) at other
than the calculated high-symmetry sites has been determined using a separate model
calculation, described in Section 3.1 and Ref. [Tom 91b], which relates End to the total
charge density of the graphite host at the Pd adsorption site. This energy functional
is assumed to be universal and reproduces End in high-symmetry sites accurately. In
this calculation, I assumed a perfectly rigid substrate and a monatomic Pd tip. A
top view of this tip at the graphite hollow site is shown schematically in the inset of

Fig. 3.6(b).

 

30-

 

 

iriio‘9 N)
N
O

 

 

 

 

 

 

 

 

 

 

Figure 3.5: (a) Calculated force F between a Pd AFM tip with one apex atom
and the surface of hexagonal graphite, as a function of the Pd-graphite distance 2.
The solid and dashed lines correspond to the sixfold hollow (H) and the on-top (T)
sites, respectively. An enlarged section of the graph near the equilibrium is shown
in the inset. A second inset shows the adsorption geometry in top view. A possible
trajectory of an AFM tip along x is shown by arrows. (b) Valence charge density
of the Pd/ graphite system. The results of the LDA calculation are for the on-top
adsorption site near the equilibrium adsorption distance 2”, and are shown in the
:2 plane perpendicular to the surface. The ratio of two consecutive charge density
contours p(n + l)/p(n) is 1.2. (From Ref. [Zho 91b]. ©Les Editions Physique)

43

 
   
      

    

III

I III].

‘ n“. 0.094,,”
\

. “3.0.21, ”‘

  
     
   
 

   
  
  

. I, I, r/
l I ’r
‘ 0 9% ',,’I/,'
o ’I I”
a 0, III

 
 
   
  

r"!
I I . . . I I!
I,’I,’;I,,, ‘. o ”obi/n.
-n\| s..: I ” “n so I II:
-\\\“| . 5’07" us‘.°.’,”I/r
\| O I .e
‘ M35 ’1" . .' '3

   

     

 

 

 

 

 

W \WW\
T 1T 11'] . T. T. H T T H
0.00 2.50 5.00 7.50 10.00 12.50

x(K)

 

Figure 3.6: (a) Surface corrugation Az experienced by a monatomic Pd AFM tip
scanning the my surface plane of rigid graphite under the applied AFM load (per
atom) fm = 10‘8 N. (b) A2 (with respect to the “H” site) along the surface .1:
direction for a “sharp” l-atom tip, for fa. = 10'” N (dotted line), 5 x 10'9 N
(dashed line), 10"8 N (solid line), and 2 x 10'8 N (dash-dotted line). The inset shows
the geometry of the tip-graphite system in top view. The AF M tip is shown above
the hollow site, and the shaded area represents the Pd atom. (From Ref. [Zho 91b].
©Les Editions Physique)

44

Fig. 3.6(b) shows the AF M corrugation A2(.r) for different loads fat. The tip
trajectory along the surface zit-direction, shown by arrows in the inset of Fig. 3.6(b),
contains the “T” and “H” sites and yields the largest corrugation. As discussed
above, the favored surface site changes with changing load. For the sake of simple
comparison, I set A2(hollow site)=0 in Fig. 3.6(b) and obtain a sign change of A2 near
for: = 2.5 x 10‘9 N. These calculations show that atomic resolution in the constant
force mode in the AFM, corresponding to A22,0.05 A, requires loads fcx¢i5 x 10‘9 N.
Since the corrugation A2 along a trajectory connecting adjacent “T” sites is very small
(see Fig. 3.6(b) ), the observation of individual carbon atoms is unlikely, which has

been confirmed by the experiment [Mey 88].

For an n-atom tip, which is commensurate with the substrate, the average load per
AFM tip atom is fat = Fat/n and the equilibrium tip position 2 can be estimated
from Eq. (3.5). It should be noted that under certain conditions, such a “dull” multi-
atom tip can still produce atomic corrugation for fat similar to a monatomic tip.
This is the case for an ideally aligned tip with a close-packed (111) surface, since the

unit cells of Pd and graphite are nearly identical in this case.

The range of applicable loads Fm is limited by the condition that substrate dis-
tortions near the AF M tip should be small and remain in the elastic region. Since full-
scale LDA calculations of local AFM-induced distortions of a semi-infinite graphite
surface are practically not feasible, I adopt the following approach. First, I use con-
tinuum elasticity theory [Lee 89], with elastic constants obtained from ab initio calcu-
lations [Tom 89], to determine the relaxation of carbon atoms at the graphite surface
in response to the AF M load applied through an AF M tip. [This continuum approach
is applicable in the linear response regime and has been successfully used previously

to calculate local rigidity, local distortions and the healing length of graphite near

45

an AFM tip and near intercalant impurities [Tom 89, Lee 89]. In a second step, the
atomic structure in the total charge density of the Pd/ graphite system is regained

from a superposition of atomic charge densities given by LDA-atom calculations.

The semi-infinite system of graphite layers is characterized by the interlayer spac-
ing d, the in-plane C-C bond length dc-c, the flexural rigidity D, the transverse rigid-
ity K (proportional to C44) and c-axis compressibility G (pr0portional to C33)[Tom 89,
Lee 89]. The LDA calculations for undistorted graphite yield d = 3.35 A and
dc—C = 1.42 A, in excellent agreement with experiment [Zab 89]. In the contin-
uum calculation, I further use D = 7589 K, K = 932 KA'2 and G = 789 KA'“
which have been obtained from calculated graphite vibration modes [Tom 89] and

the experiment [Zab 89].

The total charge density of the graphite surface, distorted by a Pd AF M tip,
is shown in Fig. 3.7. A comparison of charge density contours with results of the
self-consistent calculation in Fig. 3.5(b) proves a posteriori the applicability of the
linear superposition of atomic charge densities. Fig. 3.7(a) shows that the substrate
distortion in response to a monatomic AFM tip at a load Fm = fez: = 10’9 N is
moderate. According to Fig. 3.6(b), the corresponding corrugations during an AF M
scan with this load are A2 z 0.03 A and thus below the limit of detection. A
larger load Fm = [at = 5.0 x 10"9 N leads, according to Fig. 3.6(b), to marginally
detectable corrugations A2 a: 0.06 A, but distorts graphite much more, as shown in
Fig. 3.7(b). ‘

For larger applied forces Fm = f“. > 5.0 x 10’9 N, which would lead to sizeable
corrugations, the local distortions of graphite exceed the elastic limit. A rough es-

timate of these distortions, based on continuum elasticity theory, indicates that for

f mZS x 10‘9 N, the distance between graphite layers near the AFM tip approaches

L

46

 

 

      
 

>>X TX ’ i091 "

    

@X ‘

 

 

Figure 3.7: Total charge density p of the monatomic Pd AFM tip interacting with
the elastic surface of graphite near the hollow site. Contours of constant p are shown
in the :2 plane perpendicular to the surface, for (a) Fm = fest = l x 10'9 N and for
(b) Fm = f,“ = 5 x 10'9 N. The ratio of two consecutive charge density contours
p(n + l)/p(n) is 1.4. The location of the applied load acting on the Pd atom is
indicated by a triangle. (From Ref. [Zho 91b]. ©Les Editions Physique)

47

the value of intra—layer C-C distances. This new diamOnd-like bonding geometry leads
to a rehybridization of carbon orbitals and will necessarily result in an irreversible
substrate deformation. An independent estimate of the critical AFM force for this
plastic deformation can be obtained from a first-principles calculation [Fah 86] of the
graphite-diamond transition as a function of external pressure. These results, corre-
sponding to an “infinitely extended tip”, indicate a critical force per surface atom of
fext = 10’9 N for this transition. Due to the large flexural rigidity of graphite, this
force increases by half an order of magnitude for a one-atom tip, in agreement with

the above result.

A realistic AFM tip is more complex than the model tip discussed above and could
consist of a micro-tip of one or few atoms on top of a larger tip. A substantial por-
tion of this larger tip could, through the “cushion” of a possible contamination layer,
distribute the applied load more evenly across a large substrate area, reduce the large
curvature near the tip (see Fig. 3.7(b) ) and increase the minimum interlayer sepa-
ration. This effect would increase the upper limit of applicable loads fm compatible

with elastic substrate deformations and would lead to atomic resolution [Mey 88].

Summarizing these results, I used ab initio calculations to determine corrugations
observable in Atomic Force Microscopy of graphite. I found that in the constant-
force mode, atomic resolution is marginally possible for AFM loads (per atom) close
to 5 x 10’9 N. For smaller loads, the corrugation A2 is too small to be observed.
For loads which are too large, graphite deformations exceed the elastic limit and

subsequently result in the destruction of the substrate.

48

3.3 First-principles theory of atomic-scale fric-
tion

In this Section, I present a predictive first-principles theory of atomic-scale friction.
I use ab initio results of the interaction energy between a perfect layer of atoms and
a graphite layer, as a function of their relative distance . I use this information to
determine the friction coefficient u between these systems in contact. My results
include quantitative predictions of u as a function of the external force, a qualitative
explanation of the increase of u with increasing surface roughness, and a qualitative
explanation of the decrease of u with increasing relative velocity between the objects

in contact in case of sliding friction.

For a microscopic understanding of the friction process, I first consider the motion
of an atomic layer along the surface. This layer is in registry with the substrate and
experiences an external force per atom fez, normal to the surface. At each site, the

equilibrium adsorption length 2 is given by the condition

fest = ‘%E0d(z)- (3.6)

I consider a straight trajectory along the :1: direction in the surface connecting nearest
neighbor hollow sites which are separated by A3: (see the inset in Fig. 3.8(a)). The
position-dependent part of the potential energy V of the system has two main com-
ponents. These are variations of the adsorption bond energy and the work against
the external force fm applied to the adsorbate, due to variations of the adsorption

length. Hence,

V(Iifert) = Ead(za Z(¢afezt)) + fectz(xafext) _ %(fezt)1 (37)

49

 

 

 

 

 

 

 

 

10"N)
° .° . .
O r-b

 

 

 

' \
‘0 I, “ i, ‘t I
‘ , -0o1 -‘ I, | g, “ (‘1
u ‘s a “ ’I \ f
' \l H B .’ H B \‘
-0.3 1 . . -
0-0 2 5 5.0

. x(A)

Figure 3.8: (a) Potential energy V(x) of the Pd-graphite system as a function of
the position of the Pd layer along the surface x-direction, for external forces fm =
3 x 10’9 N (dotted line), 6 x 10'9 N (dashed line) and 9 x 10'9 N (solid line). The
inset shows the adsorption geometry and trajectory of the Pd layer in side view. (b)
Atomic-scale structure of the force along the surface f, (dashed line) and the friction
force f, = |f,| (solid line) for fm = 9 x 10‘9 N. (From Ref. [Zho 90]. ©American
Physical Society)

50

where I arbitrarily set the potential energy to zero at the hollow site by defining

I/O(fe:rt) = Ead(xHaz(fext)) + fext 2(IHff61‘t)‘ (3.8)

In Fig. 3.8(a) I show V(:c) for different external forces. I obtain the periodic potential
V(:r) at other than the calculated high-symmetry sites from a Fourier expansion over
the reciprocal lattice. I keep the lowest components of this expansion and determine
the expansion coefficients from the calculated values for the “T” and “H” sites. How-
ever, the calculated frictional coefficient depends only on extrema of the potential

V(:c), which is assumed to be at “T” and “H” sites.

From my calculation, In I find that the mechanical component dominates and is
only partly compensated by adsorption energy differences. As a result of variations

of V along 2:, there is a position-dependent force f3 along the a: direction, shown in

Fig. 3.8(b), which is given by

f1:(xa fest) = %V(xafext)- (3.9)

The maximum value of f, describes the static friction governing the onset of stick-
slip motion. The sliding friction on the other hand must be obtained as a weighted
average over this force from the energy dissipated in friction along 2:. I first con-
sider a conservative part of this process, corresponding to potential energy increase
AVmu(fm) = Vmu(f,¢¢) - Vm;n(fm) along Ax, which yields a positive value of ft.
The non-conservative part corresponds to a decrease of V(:c) along a: and a negative
value of f,. Friction losses AE, along Ax must not exceed the maximum increase in

the potential energy, hence

AE, _<_ AVnm. (3.10)

51

For very slow tracking velocities, any gain in potential energy during the dissipative
part of the friction process is efficiently transferred into surface phonons and electron-
hole pairs. Then, I can consider both sides of Eq. (3.10) to be equal, which corresponds
to f, = [ft], shown in Fig. 3.8(b). This is the first quantitative prediction of atomic-

scale structure in the friction force which has been observed recently [Mat 87].

The energy dissipated in friction along Ax can also be related to the average

friction force < f f > along the trajectory, as
AEf =< f;>A:1:. (3.11)

Using Eq. (3.10) for AEf, I obtain

1.

<f: >= Ax

awmr _ win

Applying the definition of the friction coefficient 14 = f f / fat, I find

_ <f,> _ AVma,

- 7:..— - ea.- ‘3'”)

In Fig. 3.9 I show u as a function of fat. I find a general increase of u with increasing
external force. The minimum in u(fm) near fez; = 5 x 10‘9 N is caused by the
switching of the minima in V(a:) from H to B, shown in Fig. 3.8(a). Near fez. =
5 x 10'9 N, the lowest order Fourier components of the potential V vanish. According
to the interpolation scheme discussed above, this leads to u = 0 (shown in Fig. 3.9). In
reality, a small nonzero value of u is expected due to higher-order Fourier components
of V (which are usually much smaller than those considered). From former results
presented in Section 3.2 and and Ref. [Tom 89], I conclude that if the external force
(per atom) exceeds 10‘8 N, the graphite surface is very strongly deformed [Mam 86]

and likely to be ruptured [Ski 71]. Since no plastic deformations have been observed

52

 

0.06

 

 

 

 

:5.

4.)

d 0.05

0

00-4

a 0.04

H

8

o 0.03

8

3;; 0.02

0

"-0

a 0.01

0.00 ‘ L l L
0 5 10 15 20 25 30
fu,(10'°N)

Figure 3.9: Microscopic friction coefficient )1 between a Pd AF M tip and graphite
as a function of load per tip atom fm. (From Ref. [Zho 90]. ©American Physical

Society)

53

in the AF M studies [Mat 87]. the applied forces were probably in the region fort <
10"8 N. For these values of fat, the calculated friction coefficient of the order p z 10‘2

agrees with the experiment [Ski 71, Mat 87].

In order to obtain a meaningful comparison with observable friction forces, I have
to make further assumptions about the macroscopic interface and the elastic response
of the substrate to external forces. In the simplest case, I consider an atomically
flat interface, where N atoms are in contact with the substrate, and neglect elastic
deformations. Then, the external force per atom fat is related to the total external

force F m by

1

feet = NFext- (3.14)

In Fig. 3.10 I use the calculated u(fm) to plot the total friction force F f for such a
perfectly flat interface consisting of 1500 Pd atoms. Since )1 increases with increasing

value of fat, the F f versus Fm relationship is nonlinear, which has also been observed

in the AFM experiment [Mat 87].

In the case of large external forces and an elastic substrate such as graphite, elastic
theory predicts [Lan 86] the substrate deformations to be proportional to Fez/,3. In
case of a spherical tip [Mat 87], the tip-substrate interface area and the corresponding
number of atoms in contact is proportional to F3153. Then, the force per atom fat is
proportional to Fez/,3. Hence for increasing external forces, variations of the effective
force per atom and of p are strongly reduced due to the increasing interface area.
This is illustrated by a dashed line in Fig. 3.10, where I used N = 1500 atoms for
Fm = 10'6 N. These results are in good agreement with the AFM results for a large

nonspecific tungsten tip with a radius R = 1500 A - 3000 A on graphite [Mat 87],

but show a slightly larger increase of the friction force than observed for the range of

 

 

 

 

20

 

 

N
0

  

’5
? 15 L
C i.
fl
v
A 10 f
to
CI. [ ,
V metros -. '
spear-sis;

5 l was m l
-' cousrmr . is: m.
. .~ ]

 

 

 

O 5 10 15 20
rux104N)

Figure 3.10: (a) Calculated macroscopic friction force F f as a function of the external
force Fm for a large object. The solid line describes a “flat” object, the surface of
which consists of 1500 atoms in contact with a rigid substrate. The dashed line de-
scribes friction of a large spherical tip and also considers the effect of elastic substrate
deformations on the effective contact area. The dotted line corresponds to a constant
friction coefficient )1 = 0.012. (From Ref. [Zho 90]. ©American Physical Society)
(b) Experimental measured average friction force as a funCtion of load on the tip as
it slides across the surface. (From Ref. [Mat 87]. ©American Physical Society)

55

external forces investigated.

This theory can be used also to explain the dependence of p on the roughness of
the interface and on the relative velocities of the two objects in contact. At a rough
surface, the number of atoms N in contact with the substrate is smaller than at a
flat surface which leads to an increase of fat and hence of p. Also, with increasing
relative velocity, the coupling between macroscopic and internal microscopic degrees
of freedom (phonons, electron-hole pairs) gets less efficient. Then, AE; < AVm“ in

Eq. (3.10) which causes a decrease of < f f > and u.

Summarizing the results in this Section, I determined the atomic-scale friction
associated with a layer of Pd atoms moving across a graphite substrate from ab initio
total energy calculations. I evaluated the friction energy caused by variations of the
chemical bond strength and work against an external normal load. The calculated
value of the friction coefficient is very small, in the order p z 10'2 for loads near
10‘8 N. This small value can be explained by the weak dependence of Pd-graphite
interaction on the adsorption site. I also found )1 to increase with load in agreement

with recent Atomic Force Microscopy experiments.

3.4 Ideal friction machines

In this last Section, I determine the upper limit of the average friction force < F f >,
originating from the potential energy barriers AV(Fm) during the “surface diffusion”
of Pd on graphite. There, I assumed that the energy gained by the tip atom is
completely dissipated into microscopic degrees of freedom. In the following, I will
address the mechanism which leads to the excitation of microscopic degrees of freedom
and hence to energy dissipation. I will make use of the first-principles Pd-graphite

interaction potential. Based on these data, I will determine atomic-scale modulations

fi'l

56

of the friction force F, along the horizontal trajectory of the Pd tip on graphite in the
quasistatic limit of relative velocity v —> 0. I will show that F f depends sensitively
not only on the corrugation and shape of the Pd-graphite interaction potential, but
also on the specific construction parameters of the F FM. The latter point is of utmost
importance if friction forces obtained with different microscopes are to be compared

to one another.

Two models of a Friction Force Microscope are shown in Fig. 3.11. In both models,
the suspension of the tip moves quasistatically along the surface x-direction, with its
position 2);; as the externally controlled parameter. The tip is assumed to be stiff with
respect to excursions in the surface y-direction. I restrict the discussion to the case of
tip-induced friction and assume a rigid substrate, which applies for friction measure-
ments on graphite [Mat 87, Zho 90]. In the “maximum friction microscope” [Zho 90],
the full amount of energy needed to cross the potential energy barrier AV along As
is dissipated into heat. This process and the corresponding friction force can be ob-
served in an imperfect Atomic Force Microscope which is shown in Fig. 3.11(a). A
vertical spring connects the tip and the external microscope suspension M. The hor-
izontal positions of the tip and the suspension are rigidly coupled, :c, = 2M = :c. For
0 < :1: < Ax/ 2, a constant load Fm on the tip is adjusted by moving the suspension
up or down. For Ant/2 < a: < Ax, however, the tip gets stuck at the maximum value
of 2,. At Am, the energy AV stored in the spring is abruptly and completely released

into internal degrees of freedom which appear as heat.

The calculated potential energy V(:c) during this process is shown in Fig. 3.12(a).
I considered a monatomic Pd tip sliding along a trajectory connecting adjacent hollow
and bridge sites on graphite, for a load Fm = 10'8 N. In order to predict F f, I have

used the results of my previous ab initio calculation for Pd on graphite presented in

57

 

Figure 3.11: Two models for the Friction Force Microscope (FFM). In both models,
the external suspension M is guided along the horizontal surface a: direction at a
constant velocity v = dzM/dt —v 0. The load Fm on the “sharp” tip (indicated by
V) is kept constant along the trajectory 2,(:1:;) (shown by arrows). (a) A “maximum
friction microscope”, where the tip is free to move up, but gets stuck at the maximum
2, between Az/2 and Ar. (b) A “realistic friction microscope”, where the position
of the tip 2:, and the suspension 2M may differ. In this case, the friction force F J is
related to the elongation x: - 12M of the horizontal spring from its equilibrium value.

(From Ref. [Tom 91c]. ©Les Editions Physique)

58

Section 3.1 and Ref. [Zho 90], which have been conveniently, parametrized [Tom 91b].

The force on the tip in the negative :r-direction, as defined in Fig. 3.11(a), is given by

_ 5V(xM)/82M if 0 < 2M < Ax/2

Elf“) " { 0 if Ax/2 < 2M < A2: (3'15)

and shown in Fig. 3.12(b). The nonzero value of the average friction force < F f >,
indicated by the dash-dotted line in Fig. 3.12(b), is a clear consequence of the mech-
anism which allows the tip to get stuck. F f(xM) is a non-conservative force since
it depends strongly on the scan direction. In absence of the “sticking” mechanism,
F f is given by the gradient of the potential everywhere, as indicated by the dashed
line in Fig. 3.12(b). It is independent of the scan direction and hence conserva-
tive. In this case, F ; inhibits sliding for 0 < 2354 < Az/2 and promotes sliding for

Ax/2 < 2M < Ax, so that the friction force averages to zero.

A more realistic construction of Friction Force Microscope is shown in Fig. 3.11(b).
In this microscope, the AFM-like tip-spring assembly is elastically coupled to the
suspension in the horizontal direction, so that 2:, may differ from 2M. For a given
1:5, the tip experiences a potential V(:r¢,2¢) = V,-m(.r¢,2¢) + Fm 25 consisting of the
tip-surface interaction V5,“ and the work against Fm. The tip trajectory 2¢,m,n(a:t)
during the surface scan is given by the minimum of V(:r,, 2,) with respect to 2,. For

this trajectory, V(:c¢) = V(:1:5, 2¢,m,-,,) represents an effective tip-substrate potential.

This potential V(:c¢) depends strongly on Fm and is corrugated with the period-
icity of the substrate due to variations of the chemical bond strength and of 25min.
It is reproduced in Fig. 3.13(a) in the case of a monatomic Pd tip on graphite and
Fm = 10’8 N. The corrugation of the potential V(;1:t) will elongate or compress the

horizontal spring from its equilibrium length which corresponds to 2:, = 3M. The

 

 

V(X)(eV)

 

 

 

 

 

 

 

 

 

 

 

1.6 t(b) I r 3 ,
1.0 ”- -(
A <r,>
z 05 b o-e-e -e-e-e-e-e ..... g
'1'
o 0.0 ,
V4 ‘ a
V \ I,
[gr-0.5 - x‘ .’ .
“1.0 ' ‘\‘-"I’ .1
-1.5 l n 1

0 142$ 2 A 3
1(3)

Figure 3.12: Potential energy of the tip V(:r) (a), the friction force F,(:1:) and the
average friction force < F I > (b) in the “maximum friction” microscope with a
monatomic Pd tip on graphite. The arrows indicate the tip trajectory corresponding
to a relaxed vertical position 2. for a constant load on the tip Fm = 10"8 N. (From

Ref. [Tom 91c]. ©Les Editions Physique)

60

 

 

0.8

(a)

F... 10'8 N

0.6 P

 
  

 

0.4 " d

V(xt )(eV)

 

 

 

 

 

 

 

 

/// //

 

 

 

 

 

 

it}? zAxs

MA)

Figure 3.13: Microscopic friction mechanism in the “realistic friction microscope”.
The calculations are for a monatomic Pd tip on graphite and Fm = 10"8 N. Results
for a soft spring, giving nonzero friction, are compared to a zero-friction microscope
with a hard spring. (a) Potential energy of the tip V(z.). (b) A graphical solution of
Eq. (3.19) yielding the equilibrium tip position at the intersection of the derivative of
the potential 6V(2g)/63¢ (dashed line) and the force due to the horizontal spring Ff.
Solid lines, for different values of 1:", correspond to a soft spring with c = 10.0 N/ m,
and dashed lines correspond to a hard spring with c = 40.0 N / m. (c) The calculated
equilibrium tip position z¢(zM). (d) The friction force F, as a function of the FFM
position 254 and the average friction force < F f > for the soft spring (dash-dotted
line). (From Ref. [Tom 91c]. ©Les Editions Physique)

61

friction force is given by
PAW) = -c(a=: - m). (3.16)

where c is the horizontal spring constant. The total potential energy VM of the system

consists of V(:c¢) and the energy stored in the horizontal spring,
1 2
V,O¢(:r¢, 3M) = V(a:t) + §c(a:t — 2M) . (3.17)

For a given horizontal position ml of the FFM suspension, the equilibrium position

of the tip at, is obtained by minimizing Vtot with respect to mt. This gives

aviot _ aV(I¢)

 

6:5: -— 7;:— + C(51); — xM) = 0 (3.18)
or, with Eq. (3.16),
_ _ (ii/(1%)
Ff — CED; + cxM — 61¢ . (3.19)

A graphical solution of Eq. (3.19) is shown in Fig. 3.13(b) and the resulting relation
x¢(x M) is shown in Fig. 3.13(c). If the force constant c exceeds the critical value cc,“ =
—[02V(:r¢)/8z?]m,-,,, which for Pd on graphite and Fm = 10‘8 N is cc,“ = 23.2 N /m,
I obtain a single solution x, for all 354. This situation is indicated by the dotted
line in Figs. 3.13(b) and 3.13(c) for a hard spring with c = 40.0 N /m. The friction
force F, is given by Eq. (3.16) and shown by the dotted line in Fig. 3.13(d). F; is
independent of the scan direction and hence conservative, resulting in < F f >= 0.
Hence no friction should occur in the AFM, which is the limiting case of an F FM for

c—voo.

A more interesting case arises if the horizontal spring is soft, c < cm}. This sit-

uation is shown by the solid line in Figs. 3.13(b) and 3.13(c) for c = 10.0 N/m. In

62

this case, the solution :r.(1:M) of Eq. (3.19) displays a sequence of instabilities. These
instabilities lead to a stick-slip motion of the tip with increasing mM, similar to “pluck-
ing a string”. The hysteresis in the x,(:cM) relation, shown in Fig. 3.13(c), results in
a dependence of the force F I on the scan direction. The friction force F f(:r M) in this
case is shown by the solid line in Fig. 3.13(d). It is a non-conservative/dissipative
force and averages to a non-zero value < F f >= 3.03 x 10'10 N, given by the dash-
dotted line. The energy released from the elongated spring into heat is represented

by the shaded area in Fig. 3.13(d).

The present theory predicts occurrence of friction only for very soft springs or a
strongly corrugated potential V(a:¢). The latter fact can be verified experimentally
since the corrugations AV(1:¢) increase strongly with increasing applied load [Zho 90].
Consequently, for a given c, the friction force is zero unless a minimum load Fm is
exceeded. On the other hand, for a given Fm, no friction can occur if c exceeds a

critical value can-(Fest).

A similar situation occurs during sliding between large commensurate flat sur-
faces of A on B. In that case, c is given by the elastic constants of A at the inter-
face [Toml 29], hence can not be changed independently. Since c is rather large in
many materials, zero friction should be observed for moderate applied loads in the
absence of wear and plastic deformations. For a multiatom “tip” which is commen-
surate with the substrate, the tip-substrate potential is proportional to the number
of tip atoms at the interface, n, and so is the critical value cc,“ for nonzero friction.
In this case, the effective FFM spring depends both on the external spring and the
elastic response of the tip material. The inverse value of cc,“ is given by the sum
of the inverse values of the corresponding spring constants. For a large tip which is

incommensurate with the substrate, no friction should occur [McC 89].

63

The average friction force < Ff > as a function of the load Fm and the force
constant c is shown as a contour plot in Fig. 3.14 for a monatomic or a larger com-
mensurate Pd tip on graphite. Clearly, the applicable load range is limited by the
underlying assumption of contact without wear. This figure illustrates that not only
the friction force F f, but also the friction coefficient )1 =< F f > /Fm, depend
strongly both on the interaction potential between the two materials in contact and
on the intrinsic force constant c of the Friction Force Microscope. This clearly makes
the friction force dependent on the construction parameters of the FFM. There is also
one advantage in this fact: c can be chosen in such a way that nonzero friction occurs

even at small loads Fest-

Summarizing the results in this Section, I calculated the atomic-scale modulation
of the friction force and the corresponding stick-slip motion at the interface during
the relative motion between Pd and graphite. I proposed two idealized versions of
a Friction Force Microscope. I showed that the friction force depends not only on
the Pd-graphite interaction potential, but even more critically on the construction

parameters of such a microscope.

3.5 Conclusions

In this Chapter, I have studied theory of atomic force microscopy, using Pd AFM tip

and graphite substrate as the model system.

I have calculated interaction between Pd and graphite using ab initio density
functional theory within local density approximation, with Pd atom placed above
hollow and on-top sites of graphite. The precise calculations show that Pd—graphite
interaction is characterized by strong close shell repulsion at small distances, which

strongly favors hollow site, and very weak chemical bond at larger distance with on-

 

64

c/n (Nm—l)

 

 

 

F.,./n (10‘5N)

Figure 3.14: Contour plot of the average friction force < F f > between a Pd tip
and graphite, as a function of the load Fm and the force constant c. All forces and
the force constant are normalized by the number of tip atoms 71 in contact with the
substrate. (From Ref. [Tom 91c]. ©Les Editions Physique)

65

top site slightly favored. The total energy results at high symmetry points can be
well mapped to a parameterized form, as a function of graphite host charge density.
The optimum combination of ab initio and parameterized results provides a precise

and complete picture of Pd/ Graphite interaction.

Using these results, I have studied the limit of resolution of Pd AFM tip on the
graphite surface. The calculated trajectories of Pd probing graphite surface under
different loads reveals that, the minimum load on Pd is about 5 x 10’9N to produce
observable corrugation. Using continuum elasticity theory, I have estimated that the
load exceeding 10’3N leads to too much deformation and is thus destructive to the
surface. So the atomic resolution of graphite from Pd AFM tip is only marginally

achievable.

The atomic scale friction between an Pd tip and graphite surface is obtained
from the variation of Pd tip energy along the scanning trajectory. There are two
contributions to this energy, one is from the work against the load when corrugation
is none zero, another comes from adsorption energies difference at different sites.
Assuming all the kinetic energy gained by the AFM tip dissipates into heat, I obtain
the first quantitative prediction of friction force using ab initio techniques. The
calculated friction coefficient is in the order of 10'2 for loads near 10’8 N, and it
increases with increasing loads. The calculation is in good agreement with recent

experiment.

I further investigate the possible mechanisms of energy dissipation during the fric-
tion process. I proposed two models for friction force microscope (FFM). In the more
realistic version, I was able to calculate the atomic scale modulation of the friction
force, and thus explain the stick-slip motion of the tip. I found the averaged friction

force depends not only on the tip-surface interaction potential, but also strongly on

66

the construction parameters of the FFM.

Chapter 4

Structure and Dynamics for H

Loaded Pd

As I already mentioned in Chapter 1, the structure and dynamics of H loaded Pd is of
great importance to both basic science and technology. While the ultimate goal is to
perform the corresponding study completely from first principles, the computational
involvement is presently prohibitive for such an endeavor. Therefore, I base the study
of the Pd—H system on the Many-Body Alloy (MBA) Hamiltonian, which was derived
in Chapter 2, with parameters reproducing the results of an LDA calculation for this
system [Tom 86b, Sun 89, Tom 91a]. In Section 4.1, I use the Hamiltonian to calculate
the structural properties of bulk Pd and PdH, as well as clean and hydrogen covered
Pd surfaces. The results of the surface studies, which have no additional adjustable
parameters, are found to be in good agreement with corresponding LDA calculations

and experimental data.

In Section 4.2, I use lattice dynamics to calculate the dispersion relations for
bulk Pd and PdH, and compare the results to experimental data. The effect of H
adsorbates on the phonon spectra of Pd (001) and Pd (110) surfaces is investigated

by performing the corresponding lattice dynamics calculations.

In Section 4.3, dynamical properties of Pd-H systems at finite temperature are

67

!.

68

investigated using molecular dynamics. I present results for the equilibrium struc-
tural and elastic properties of Pd at different temperatures and hydrogen concentra-
tions, and show quantitative results for the mechanical failure and crack formation in
hydrogen-free and hydrogen-loaded Pd due to a uniaxially applied load. I will show
that a careful analysis of these results can provide valuable insight into the mechanism

of “hydrogen embrittlement”.

This Chapter contains material which has appeared in the following three publi-

cations:

e W. Zhong, Y. S. Li, and D. Tomanek, Effect of Adsorbates on Surface Phonon
Modes: H on Pd(001) and Pd(110), Phys. Rev. B44, 13053 (1991).

e W. Zhong, Y. Cai, and D. Tomanek, Mechanical Stability of Pd-H Systems: A
Molecular Dynamics Study, Phys. Rev. B 46 8099 (1992).

e W. Zhong, Y. Cai, and D. Tomanek, Computer Simulation of Hydrogen Em-
brittlement in Metals, Nature 362 435 (1993).

4.1 Application of the Many-Body Alloy Hamil-
tonian to the Pd—H system

In this Section, I determine the parameters of the MBA Hamiltonian for the Pd—H
system, based on previously published ab initio results[Tom 86b, Sun 89, Tom 91a]
for the equilibrium structure and cohesive energy of bulk Pd and PdH. It should be
noted that this procedure does not introduce any free parameters. In the following,
I will apply this Hamiltonian next to determine structural properties and vibrational

spectra of clean and hydrogen covered Pd surfaces.

 

69

Table 4.1: Interaction parameters used in the Many-Body Alloy Hamiltonian for the
Pd-H system.

 

Interaction qag pag rag,o(A) 653(6V) {05(eV)
Pd-Pd 3.40 14.8 2.758 0.08376 1.2630
H—H 3.22 5.28 2.300 0.1601 0.9093
H-Pd 2.20 5.50 1.769 0.6794 2.5831

 

 

4.1.1 Construction of the Many-Body Alloy Hamiltonian

In the MBA Hamiltonian, each of the H—H, H—Pd and Pd—Pd interactions is char-
acterized by a set of five parameters: {0,63,24,12 and r0 (four of these parameters
are independent). The Pd—Pd interaction is obtained from the previous ab initio
calculation[Tom 91a] of the cohesive energy Eco}, as a function of the lattice constant
a for bulk Pd. I consider nearest neighbor interactions only and obtain a simplified

expression for the bulk cohesive energy,

 

7‘?de _1)]}1/2

r0,Pde
rPde _1)]

r0,Pde

E...(Pd bulk) = —{ 2.... 53.....epr—2qp.,p.(

 

+ 21>qu ‘53.?de eXPl-P( (4-1)

Here, Zbuu‘ = 12 for the fcc structure and mpdpd = a\/2/2 is the nearest neighbor
distance. The calculated cohesive energy Eco}, of bulk Pd as a function of the lattice
constant a is given by the solid line in Fig. 4.1 and compared to corresponding LDA

results of Ref. [Tom 91a]. The parameters used in Eq. (4.1) are given in Table 4.1.

The parameters for H—H interaction can be determined in a similar way as those
for the Pd—Pd interaction, by mapping the MBA Hamiltonian to the ab initio results

for molecular hydrogen. The corresponding parameters are given in Table 4.1.

To determine the parameters for the Pd—H interaction, I apply the MBA Hamil-
tonian to bulk PdH. Considering the nearest-neighbor Pd—H, H—H and Pd-Pd in-

teractions, I formally decompose the cohesive energy of the PdH crystal with N aCl

 

70

 

0.5 '-

0.4 '-

 

0.3

AEcoh (eV)

0.2

0.1

 

 

 

 

 

0.0
3.7 3.8 3.9 4.0 4.1 4.2 4.3 4.4 4.5

a (A)

Figure 4.1: Cohesive energy changes AEcoh = E60,, - Ecohp in bulk Pd and PdH as a
function of the lattice constant a. Values obtained using the MBA Hamiltonian for
Pd (solid line) and PdH (dashed line) are compared to LDA results of Ref. [Tom 91a]

for the corresponding systems, given by e and 0. (From Ref. [Zho 91a]. ©American
Physical Society)

 

71

structure as

Ecoh = Ecoh(H) "l" Ecoh(Pd)- (4.2)

The binding energies of H and Pd atoms in this structure are given by

 

Ecoh(H) = -{ ZH(Pd)§<i.PdH exPl-2QPJ,H( rm” " 1)]

7‘0,PdH

 

1/2
+ZH( (1)5333 CXPI-2QH.H( rm” - 1M) (4-3)

V

THE
7‘0,PdH

 

+{ ZH(Pd)€3,de expl‘PPd,H( - 1)]

 

+ZH(H)€ciiHu expl—PH.H( rm! - 1”)

 

 

 

 

T0,}!!!
and
7'
E...(Pd) = —{ ZP.(H)53,p.Hexpl-2qp.,n( Pd” -1)]
7‘0.PdH
2 7‘?de ”2
”pawn...” expl—2qp.,p.( - In} (4.4)
70,13de
,.
+{ ZPd(H)€gPdHexPl"PPd.H( Pd” *1”
7‘0,PdH
,.
mama...epr-pp.,p.( ”P" —1)1}.
r0,Pde

Here, Zpd(H) = 6, Zpd(Pd) = 12 are the respective numbers of H and Pd nearest
neighbors of a Pd atom, and ZH(Pd) = 6, ZH(H) = 12 are the numbers of H
and Pd nearest neighbors of a H atom. 1'0,de = a/2, rapdpd = a\/2/2 are the
Pd — H and Pd — Pd nearest neighbor distances and a is the lattice constant. The
values of (0.13.111, pde-I, qufl, and 6,1,1,de have been determined by reproducing LDA
results of Ref. [Tom 91a] for Each of bulk PdH and are given in Table 4.1. The

corresponding results obtained with the MBA Hamiltonian are given by the dashed

 

72

line in Fig. 4.1, together with the LDA results given by the data points. The good
agreement between results for the bulk systems based on the MBA method and LDA
calculations indicates that the model Hamiltonian has sufficient flexibility to describe

energy changes accurately.

4.1.2 Structure and stability of clean and hydrogen covered
Pd(001) and Pd(110) surfaces

In the previous Subsection, I showed that an appropriately parametrized MBA Hamil-
tonian reproduces the bulk equilibrium properties very accurately. In this Subsection,
I will test the applicability of the MBA Hamiltonian to surfaces, specifically to the
calculation of surface energies, surface relaxations, and adsorption energies of H on
(001) and (110) surfaces of Pd. These results will be compared to experimental data

and to ab initio data of Refs. [Tom 86b] and [Tom 91a].

The surface energy of clean Pd surface can be obtained from the cohesive energy

of the bulk and that of an n-layer slab using[Tom 91a]
E', = %[ECOh(Pd slab) — nEcoh(Pd bulk)]. (4.5)

This expression is relating values of the slab energy Ecoh(Pd slab) and the correspond-
ing surface energy E, per atom in a surface or a layer. Using this equation together
with MBA Hamiltonian, one can easily determine the surface energy for a given Pd

surface. Multilayer surface relaxations can then be determined by minimizing E,.

I have calculated the surface energy for Pd slabs with both (001) and (110) sur-
faces. MBA results for the change of the surface energy AE, due to a relaxation Adm
of the topmost interlayer distance at the (001) and (110) surfaces of Pd are given by

the solid lines in Figs. 4.2(c) and (d), respectively. For the sake of simple comparison

with the LDA results of Ref. [Tom 91a] [given by the data points in Fig. 4.2(d)], the

73

MBA calculations have also been performed for a 3-layer Pd slab. My results indicate
a surface contraction which increases with a decreasing coordination number of the

surface atoms, in good quantitative agreement with the LDA data.

It is easy to handle very thick slabs with the MBA Hamiltonian, and I will con-
centrate on much thicker slabs with n = 25 in the following. This slab thickness is
more than sufficient to guarantee that the two slab surfaces do not interact and that

the atoms in the middle of the slab are truly in a bulk environment.

These numerical results indicate that the surface energy of Pd(110) (0.73 eV / atom)
is much higher than that of Pd(001) (0.48 eV/atom). As discussed in more detail
in Ref. [Tom 91a], this can be tracked back to the narrowing of the effective Pd 4d
bandwidth, or an increasing number of “dangling bonds”, with decreasing coordi-
nation number. Since precise experimental data for these surface energies are not
available, I will only compare these results to recent LDA calculations[Tom 91a]. For
Pd(001), the surface energy value E, = 0.48 eV/atom is very close to the the LDA
result[Tom 91a] of 0.49 eV/atom. The MBA value E, = 0.73 eV/atom for a Pd(110)
surface is significantly lower than the reported LDA result of 1.80 eV/ atom. The
large discrepancy between these latter results is possibly due to a less adequate basis

set and the neglect of surface relaxations in the LDA calculation.

My results for the equilibrium structure of Pd (001) and Pd (110) are summarized
in Table 4.2. The values in columns 2—4 of Table 4.2 indicate a damped oscillatory
behavior for the surface relaxations. These oscillations occur as a general phenomenon
which has been observed[Dav 83, Bar 85, Sko 87] and calculated[Tom 86a, Tom 85b,
Bob 90] in many systems. I also find the surface relaxations to be more pronounced
on the more open (110) surface than on the close packed (001) surface, in agreement

with the general observation that increasing relaxations correspond to larger surface

74

 

 

 

AEs(oV)

 

 

 

 

 

 

 

 

 

 

 

 

”0.3 1 l L n 1 m _0.3 A L n 7 n 1 n
-20 -15 -10 -5 O 5 10 15 20 '20 -15 -10 -5 0 5 10 15 20
Ad“ (2) Mrs (K)

Figure 4.2: Schematic side View of the Pd(001) (a) and Pd(110) (b) surfaces showing
the definitions of the hydrogen adsorption height h and the interlayer spacings du, d”
and J34. Surface energy changes AB. = E, - Em for clean and H-covered Pd(001)
(c) and Pd(110) ((1) surfaces, as a function of the first interlayer spacing du. Values
obtained using the MBA Hamiltonian for clean (solid line) and H—covered (dashed
line) surfaces are compared to the LDA results of Ref. [Tom 91a] for the corresponding
systems, given by e and 0. (From Ref. [Zho 91a]. ©American Physical Society)

75

Table 4.2: Relaxations at clean and hydrogen covered Pd surfaces.

 

 

 

Surface Ad12(%) Ad23(%) Ad34(%) Reference
Pd(001) -2.2 0.2 0.0 Present Work
Pd(110) —5.2 0.7 —0.3 Present Work
Pd(001) + p(l x 1) H 4.0 0.7 0.0 Present Work
Pd(110) + p(1 x 1) H 1.3 1.3 0.0 Present Work
Pd(110) —6.0:l:2 1.0:t2 — [Bar 85]
Pd(110) —5.1:l:1.5 2.9i1.5 — [Sko 87]
Pd(110) + (2 x 1) H —2.2:l:1.5 2.9:l:l.5 — [Sko 87]

 

energies. As shown in Table 4.2, the surface relaxations, calculated using the MBA
technique, are in gratifying agreement with both experimental data[Bar 85, Sko 87]
and LDA results[Tom 91a, Bob 90]. This good agreement confirms that the MBA
Hamiltonian can accurately describe the equilibrium structure and energy changes of

clean Pd(001) and Pd(110) surfaces.

Finally, I use the MBA Hamiltonian to determine the binding energy of hydrogen
at different adsorption sites on the Pd(001) and Pd(110) surfaces. These calculations
yield the preferential adsorption site and adsorption height, both of which are ac-
cessible to experimental verification. In order to simplify the comparison with the
LDA calculations of Refs. [Tom 86b] and [Tom 91a], I performed all calculations for a
3—layer Pd slab which was covered by H on both sides. I found the adsorption energy
of H on Pd to be typically of the order E“; z —3 eV. Since the binding energy of
a H; molecule is only De = 4.75 eV, the dissociation probability of a H; molecule

approaching the Pd surface is high. For this reason, the present investigations has

been limited to atomic hydrogen on Pd(001) and Pd(110).

The equilibrium structure of hydrogen covered surfaces can be determined in the
same way as that of the clean surfaces described above. My results for the surface

energies and multilayer relaxations of hydrogen covered surfaces are summarized in

76

Fig. 4.2 and Table 4.2. Hydrogen atoms are assumed to occupy the equilibrium
sites during the surface relaxation. The calculated changes of the surface energy
at hydrogen covered Pd(001) and Pd(110) surfaces are given by the dashed lines
in Figs. 4.2(c) and ((1), respectively. For the sake of simple comparison with the
LDA results of Ref. [Tom 91a] [given by the data points in Fig. 4.2(d)], also these
model calculations have been performed for a 3~layer Pd slab. My results indicate an
expansion of the hydrogen covered surfaces, in agreement with the LDA data. The
reversal of surface contraction obtained for clean surfaces can be explained by the

saturation of Pd dangling bonds by H atoms.

In Fig. 4.3, I display the adsorption energy of H on Pd(001) as a function of the
adsorption height h for the on-top site, the bridge site and the hollow site. Results
obtained using the MBA Hamiltonian, shown by lines, are compared to LDA data of
Ref. [Tom 86b], given by the data points. As discussed above, I represent the substrate
by the same 3-layer Pd slab as used in the LDA calculation. Since it is only energy
differences which are relevant for the preferential adsorption sites and the vibration
frequencies, I concentrate on adsorption energy changes AB“ = End - Emu with
respect to the adsorption energy of H at the equilibrium site. The corresponding
results are shown in Fig. 4.3(b). I find the fourfold hollow site to be the equilibrium
site on Pd(001), in agreement with the experimental data of Ref. [Beh 80] and the
LDA results of Ref. [Tom 86b]. Also the very small equilibrium adsorption height ho =
0.18 A above the first Pd layer is in good agreement with the LDA value[Tom 86b]
of 0.24 A and the experimental result[Bes 87] of 0.30 A. These small adsorption
heights result from the very small atomic radius of hydrogen. More important, the
MBA Hamiltonian gives H-Pd interaction potentials (and consequently vibration

frequencies) in close agreement with the LDA data.

 

 

 

 

A 1.5 - 0', . -
e O .‘
1 o ' “ q...‘ cl
8 ‘ I
U ‘s I
a . ‘Q.’

 

 

 

'-o.5 0.0 0.5 1.0 1.5 2.0

h (A)

Figure 4.3: (a) Schematic top view of the Pd(001) surface and the assignment of ad-
sorption sites: hollow (x), bridge (0) and on—top(D). (b) Adsorption energy changes
AB“ = E“ - E.“ (with respect to the equilibrium adsorption energy) as a function
of the hydrogen adsorption height h on Pd(001). Results based on the MBA Hamil-
tonian for the hollow site (solid line), the bridge site (dashed line) and the on—top site
(dotted line) are compared to LDA results of Refs. [Tom 86b] and [Tom 91a], given
by the data points. (From Ref. [Zho 91a]. ©American Physical Society)

78

 

 

 

 

 

Figure 4.4: (a) Schematic top view of the Pd(110) surface and the assignment of
adsorption sites: hollow (x), long-bridge (o), short—bridge (o) and on-top (D). (b)
Adsorption energy changes AB“ = E“, — 5.1.1.0 (with respect to the equilibrium
adsorption energy) as a function of the hydrogen adsorption height h on Pd( 110).
Results based on the MBA Hamiltonian for the hollow site (solid line), the long—
bridge and the short-bridge sites (dashed line) and the on-top site (dotted line)
are compared to LDA results of Ref. [Tom 91a], given by the data points. (From
Ref. [Zho 91a]. ©American Physical Society)

 

79

Adsorption energies for H on Pd(110) are shown in Fig. 4.4 for the on-top, the
short-bridge, the long-bridge and the hollow site. In analogy to Fig. 4.3, the continu-
ous lines represent the MBA Hamiltonian results, and the discrete data points show
the LDA values of Ref. [Tom 91a]. Also in this case, the LDA and MBA Hamiltonian
calculations have been performed for a 3—layer Pd slab. Similar to the results for
the (001) surface, the MBA potentials represent energy differences and adsorption

potentials, which are in remarkably good agreement with the LDA data.

Both the MBA technique and the LDA predict the long-bridge site to be the most
favored among the adsorption sites considered here. My calculation indicates that
in this site, the equilibrium hydrogen adsorption height is only ho = 0.09 A above
the topmost Pd layer. Recent low-energy electron diffraction (LEED)[Sko 87] and He
scattering[Rie 83, Bas 89] experiments suggest that the preferential H adsorption site
is the threefold coordinated site in the troughs on the Pd(110) surface. In this site, the
separation between the hydrogen atom and the two nearest neighbors in the topmost
Pd layer, as well as the closest Pd atom in the second layer, is 2 A. The mutually
repulsive interaction between nearest neighbor hydrogen atoms stabilizes a zig-zag
adsorption pattern in the troughs on the surface[Rie 83]. Results based on the MBA
Hamiltonian indicate that this adsorption geometry is slightly energetically disfavored
when compared to the preferential long-bridge site. While the MBA Hamiltonian
clearly can not resolve such minute energy differences, I find that the predicted H—Pd
interaction potentials are in remarkably good agreement with ab initio calculations

and experimental data in view of the simplicity of the approach.

80

4.2 Phonon structure in the bulk and at the sur-
face of Pd-H systems

In this Section, I present phonon calculation results for bulk Pd and PdH, as well as
the clean and H covered (001) and (110) surfaces of Pd. As outlined in Section 2.3, the
phonon dispersion relations can be determined directly by diagonalizing the dynami-
cal matrix Dag,,w(lc.), which is basically a Fourier transformed force-constant matrix.
According to Eq. (2.26), the force-constant matrix can be determined numerically by
calculating total energy differences with respect to the atomic displacements. This
can be done even more efficiently using the MBA Hamiltonian which, as shown above,
gives reliable potential energies near the equilibrium structure. It is worthwhile to
note that the force-constant matrix accounts for effective second- and third-neighbor
interactions. This is a consequence of the many-body nature of the MBA Hamil-
tonian which correctly describes the indirect interaction between two atoms and its

mediation through a third atom neighboring the two sites.

In Fig. 4.5(a), I compare the phonon calculations for bulk fcc Pd to measured
phonon spectra along the high symmetry lines. The experimental data[Mil 71] have
been obtained using inelastic neutron scattering. As shown in this figure, the cal-
culation is in quite good agreement with the experiment throughout the Brillouin
zone. While my Hamiltonian is based on LDA calculations for static properties of
Pd, and does not contain any adjustable parameters, the good agreement with the
experimental data indicates that the MBA Hamiltonian also correctly describes the

dynamical properties of this system.

Results of a similar calculation for bulk PdH are shown in Fig. 4.5(b). The net
effect of hydrogen on the phonon spectrum is a softening on the acoustic branches.

This is partly caused by a 6% increase of the lattice constant upon hydrogen uptake

81

 

 

8.00

6.00 .

4.00 '

V0712)

2.00

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0,00.
1‘ X'XK 1" LP X'XK 1‘ L

Figure 4.5: Phonon dispersion relations for bulk Pd (a) and PdH (b). Frequencies
calculated using the MBA Hamiltonian are given by the solid lines. The experimental
data points of Ref. [Mil 71] for bulk Pd are given by (+). (From Ref. [Zho 91a].

©American Physical Society)

82

from 3.89 A in Pd to 4.12 A in PdH, as shown in Fig. 4.1. Based on the analysis
of these results, I find that the presence of hydrogen has a strong effect on the force
constant matrix <I>. The force constants describing the restoring forces acting on a Pd
atom displaced along a high-symmetry direction are reduced from c = 10.6 eV/A2 in
bulk Pd to only 6.8 eV/A2 in PdH. The effect of the the presence of hydrogen on the
bulk modulus is comparably small. I obtain B = 2.03 x 1012 dyn/cm2 for bulk Pd
which compares well with the experimental value[Kit 86] B = 1.81 x 1012 dyn/cm2
and the LDA value of 2.15 x 1012 dyn/cm“. The corresponding value in. PdH is only
3% smaller, B = 1.98 x 1012 dyn/cmz, which is in fair agreement with the LDA result

1.95 x 1012 dyn/cm2 of Ref. [Tom 91a].

Since my main interest is the effect of hydrogen on the Pd modes, I do not include
the hydrogen-derived optical modes in Fig. 4.5. These modes have a very high fre-
quency and are well separated from the acoustic Pd-derived modes due to the mass
disparity of these atoms. While the hydrogen atoms can be basically thought of as
Einstein oscillators, the weak coupling between hydrogen atoms broadens the opti-
cal states to a z 3 T Hz broad band. Since the H—H nearest neighbor distance in
PdH 113,1"! = 2.06 A is much larger than the H-H interaction range (is 1 A), I have
neglected the direct H-H interaction in this calculation, but have accounted for the
dominant indirect interaction mediated by Pd atoms. A closer analysis of the opti-
cal bands reveals that two out of the three bands show no dispersion and represent
Einstein modes of hydrogen atoms with no direct coupling. The above mentioned
dispersion of the third branch reflects the degree of Pd mediated indirect interaction
between H atoms, which is almost equally strong for the first and second neighbors.

These effects in the optical band can also be seen in the measured phonon dispersion

relations of the related systems PdDo,33[Row 74] and PdTo_7[Row 86].

83

The calculated phonon spectra for the Pd(100) and Pd(110) surfaces are shown
in Figs. 4.6 and 4.7, respectively. These surfaces are represented by relaxed 25-
layer Pd slabs. Phonon spectra of the clean (001) and (110) Pd surfaces, shown
in (a), are compared to results for hydrogen covered surfaces, shown in (b). In
order to distinguish surface from bulk states, I are showing phonon bands of bulk
Pd in Figs. 4.6(c) and 4.7(c). For the sake of simple comparison, the bulk bands are
projected onto the same two-dimensional Brillouin zones in (c) as used in (a) and
(b). Note that this procedure yields continuous bulk bands at each iii-point in the
Brillouin zone. The apparent discreteness in Figs. 4.6(c) and 4.7(c) results from a
finite Iii-point sampling in a direction perpendicular to the Brillouin zone. For both
surfaces, my calculations indicate the presence of surface modes which appear either

in the bulk band gaps or are split off from the bulk band edges.

Based on the comparison of phonon spectra for Pd(001) in Fig. 4.6(a) and the
projected bulk spectra in Fig. 4.6(c), the presence of the surface introduces a soft
Rayleigh mode 51 which is substantially softer than any bulk mode in the Brillouin
zone. At the M point, this mode corresponds to vibrations of surface atoms per-
pendicular to the surface. The origin of the mode softening is a decreased interlayer
interaction at the surface, which is only partly compensated by the surface contrac-
tion. The analysis of my results indicates that the zone-edge frequency of the 51
mode increases from 3.72 THz for the unrelaxed surface to 4.11 THz for the relaxed
surface. A second surface mode S.., which is barely split from the bulk band at M,
is a transversal mode corresponding to in-plane vibrations. In addition to these soft
modes, the surface introduces a phonon mode 56 with u z 6 THz in theigap of the
bulk spectrum near X. This mode corresponds to in-plane vibrations of topmost

layer atoms, coupled to out-of-plane vibrations of second layer atoms.

 

84

 

    
        
      

6.00 ~‘~\§

A ‘ \ 5%)
N 4 00 " 3’31}:
E ‘ 2‘. l:
3’ 7&3:

       

/
ll

2.00 . ’

        

 

 
 

., J(

'8

 

I
'1
n .
>0
~11

 

v (T112)

 

 

fl} ’
7////
W

/

0']. j
,’ ///;-,.

,.,
,/tl

v (1111)

 

 

Figure 4.6: (a) Calculated phonon dispersion relations for a clean 25-layer Pd slab
with a (001) surface. (b) Corresponding results for a H-covered Pd slab (mono-
layer coverage, hollow site). (c) Bulk phonon dispersion relations of Fig. 4.5(a),
projected onto the two-dimensional surface Brillouin zone used in (a) and (b). (From
Ref. [Zho 91a]. ©American Physical Society)

 

85

 

 
     

    

     

   

/”*\

A:

      
  

\
1“
\
)I
A

h
v

Ll

’f
1))

will;
\‘i
(in.

    
 
  
  

I I
1.!

§
’ I
fl
\

\

l
l
.\
...

“L
\\\

$\

1

l l

/f.

. I")

I //

“(fitfll

\

     

\1
. .O‘A‘lu
ll"
/
”I
I
(
i

A

   

      
  

    

    

   

v('l'Hz)

   

 
 

2.00 '1

  

 

 

 

“Hill

I
\

J/

Y a

\Q

. \\ ‘\

 

 

 

141112)

 

 

Figure 4.7: (a) Calculated phonon dispersion relations for a clean 25—layer Pd slab
with a (110) surface. The experimental data of Ref. [Lah 87] are given by e. (b) Cor-
responding results for a H-covered Pd slab (monolayer coverage, long-bridge site). (c)
Bulk phonon dispersion relations of Fig. 4.5(a), projected onto the two—dimensional
surface Brillouin zone used in (a) and (b). (From Ref. [Zho 91a]. ©American Physical

Society)

 

86

The phonon spectrum of a hydrogen-covered Pd (001) surface has been shown
in Fig. 4.6(b). I have assumed a monolayer coverage corresponding to occupying
all hollow sites by H atoms. The Pd surface has again been represented by a 25-
layer slab. Similar to the bulk Pd-H system, I am mainly interested in the effect of
hydrogen on the Pd surface modes, and do not show the H-derived high-frequency
optical modes which are well separated from the Pd modes. The most striking change
in comparison to the H-free surface shown in (a) is a massive softening of the Rayleigh
mode. At M, the frequency of the 51 mode decreases from 4.11 THz to 2.77 THz
due to hydrogen adsorption. On the other hand, hydrogen does harden other modes,
such as the bulk band gap mode at u z 6 THz which has been discussed above. A
detailed analysis of the eigenstates shows that surface phonon modes with a vibration
amplitude perpendicular to the surface experience a large amount of softening. Other
surface modes with amplitudes restricted to the surface layer, such as the Se mode,
experience a hardening due to the restricted movement in presence of hydrogen atoms.
The eigenvector analysis also indicates that hydrogen adsorption also increases the
confinement of surface modes to the few topmost layers and shortens the penetration

depth into the bulk.

As done above for Pd (001), in Fig. 4.7, I compare the phonon spectra of clean and
hydrogen-covered Pd(110) surfaces to bulk phonon spectra. For the hydrogen covered
surface, I have assumed a monolayer coverage corresponding to the occupation of all
long-bridge sites (the equilibrium adsorption site) and represented the Pd surface by
a 25-layer slab.

A comparison between the calculated phonon spectra for clean Pd(110) and in-
elastic He scattering data of Ref. [Lah 87] is shown in Fig. 4.7(a). The softest surface

mode is the Rayleigh mode SI. My calculation reproduces the general features of this

87

mode quite well, but the calculated frequencies are z 10% higher than the observed
data. While I am aware that the parametrized MBA Hamiltonian may not describe
the detailed changes of metal bonding at surfaces to a very high accuracy, I can not
exclude the possibility of a slight hydrogen contamination of the Pd sample used in
Ref. [Lab 87], which is very hard to detect and would also soften the surface Rayleigh

mode.

Since Pd(110) shows the largest surface contraction, the softest surface phonon
modes are quite close to the lowest lying bulk bands throughout the surface Brillouin
zone. The strongest softening can be observed at the I7 point, where I obtain three
surface modes 51, S; and E well below the bulk band. My analysis of the eigenstates
indicates that these modes correspond to in—plane (along the surface a: and y direc-
tions) and out-of-plane (along the z direction) vibrations of .the topmost layer. The
lowest mode is an in-plane mode with an amplitude along the y direction. Similar
to the (001) surface, the surface contraction generally hardens the surface phonon
modes. This calculation shows that at the I7 point, because of the relaxation, the
lowest surface modes with topmost layer amplitudes along the a: , y, and z direc-
tions are shifted from 2.36 THz, 2.29 THz and 2.77 THz to 2.71 THz, 2.23 THz and

2.94 THz, respectively.

From the comparison between Figs. 4.7(a) and (c), one can see that the presence
of the (110) surface introduces several other vibration modes beyond the Rayleigh
mode. The bulk phonon spectrum, shown in Fig. 4.7(c), contains gaps near X, at
u z 5.5 T112, and near 17, at V a: 4.5 THz. On Pd( 110), three surface states appear
in both gaps. The 5.4 THz gap mode at X corresponds to topmost layer atoms
vibrating along the surface a: direction, and the 4.7 THz gap mode at 1" corresponds

to an analogous vibration along the surface y direction. The surface also introduces

88

high frequency states slightly above the top of the bulk bands, well visible at the
X and 17 points. These modes have considerable amplitudes deep into the bulk,
and correspond to alternating in-plane and out-of-plane vibrations on the individual

layers.

Calculated surface phonon dispersion relations for the H covered Pd(110) surface
are shown in Fig. 4.7(b). As for bulk PdH and H/Pd(001), I do not show the H-
derived high-frequency optical modes for H/Pd(110) in Fig. 4.7(c). The general trends
discussed above for the effect of hydrogen on the surface modes of the Pd(001) surface
hold also for the (110) surface. Specifically, One can observe a softening of surface
modes involving an out-of-plane motion of topmost layer atoms, and a hardening of
modes involving an in-plane vibration of surface atoms along the a: direction at X

and the y direction at 17. Other vibration modes experience very little change.

Hydrogen-induced bonding changes at the Pd(110) surface can be understood by
investigating the three lowest phonon modes SI, S; and E at the I7 point. The
51 mode, which involves in-plane vibrations of the topmost layer along the surface
y direction and out-of-plane vibrations of the second layer, has been pushed up in
frequency from 2.23 THz for the clean surface to 2.71 THz for the H covered surface.
The main reason for this hardening is a restricted freedom of motion of topmost layer
Pd atoms due to hydrogen atoms adsorbed in the long-bridge site. The S; mode,
corresponding to a vibration of topmost layer Pd atoms along the close-packed surface
:1: direction, has been softened by the presence of hydrogen from 2.71 THz for the clean
surface to 2.46 THz for the H covered surface. Finally, the E' mode, corresponding to
out-of-plane vibrations of the topmost layer, experiences the strongest softening from

2.94 THz for the clean surface to 2.56 THz for the H covered surface.

As compared to the (001) surface, phonon modes at the Pd(110) surface are less

 

89

affected by hydrogen adsorption. This is especially truefor the softening of the
Rayleigh mode. This last conclusion may not hold in case of an adsorption in the
quasi-threefold site discussed in the previous Section. In that case, I would expect
an increasing Rayleigh mode softening with increasing hydrogen coverage, which has
been observed at the F point on the related Ni(110) surface in high—resolution EELS

experiments[Leh 87].

From this calculation, one can see that hydrogen adsorption has a profound effect
on both the equilibrium structure and dynamical properties of Pd surfaces. These two
effects are closely related. The presence of hydrogen adsorbates reverses the topmost
layer contraction, and the large anharmonicity of the inter-layer interactions causes
the effective inter-layer force constant to decrease with increasing inter-layer distance.
This results in a softening of the Rayleigh surface phonon mode with an out-of~plane
amplitude in the topmost layer. Since the hydrogen-induced expansion is larger on
the more densely packed (001) surface than on the (110) surface, the Rayleigh mode

softening is also stronger on the Pd(001) surface.

At both Pd(001) and Pd( 110) surfaces, the equilibrium adsorption height of hy-
drogen is close to zero, i.e. hydrogen atoms are buried inside the topmost Pd layer,
where they affect the effective Pd-Pd force constants most. There is very little sur-
face stress associated with this adsorption site, due to the small size of the H atoms.
Surface phonon modes with in-plane vibration amplitudes are nearly independent of
surface relaxations, but are affected by the presence of hydrogen. This latter effect
results from the hydrogen-induced change of Pd—Pd force constants, and also the
restricted freedom of motion due to the extra adsorbed atoms which hardens some of

these modes.

 

90

4.3 Mechanical stability of Pd-H systems

4.3.1 Thermal expansion and the melting transition in Pd

The first nontrivial application of the present formalism is to determine the equilib-
rium volume of Pd as a function of temperature, and to study the melting transition.
In the corresponding molecular dynamics simulation, I consider an originally cubic
volume (or simulation box) containing 500 Pd atoms on an fcc lattice. Periodic
boundary conditions are used to eliminate surface effects, and the external pressure
is set to zero. I start the simulation by first equilibrating the system for 10,000 time
steps (corresponding to 2 x 10'113) in a VTN canonical ensemble which is in con-
tact with a heat bath at T = 300 K. This simulation is followed by another 20,000
time steps in the T tN canonical ensemble at zero external pressure. The equilibrium
volume V per atom is related to the volume of the simulation box, and is obtained
from a running average during the simulation. The statistical error based on the last
10,000 time steps is found to be negligible. After V has been determined for the given
temperature, the heat bath is gradually warmed up to another temperature, and the
equilibrium volume per atom is obtained from a statistical average over 10,000 time

steps at the new temperature.

My results, presented in Fig. 4.8(a), indicate a thermal expansion of the lattice
in the whole temperature range studied. In the temperature range between T =
2000 - 2050 K, the equilibrium volume per Pd atom shows a discontinuous increase,
indicative of a first-order phase transition. This phase transition corresponds to the
melting point, as can be verified by inspecting the trajectories of individual Pd atoms
at temperatures below and above the critical temperature, (shown in Figs. 4.9(b) and

4.9(c). Near the melting point, my calculations show a narrow hysteresis describing

 

91

 

N
....

N
O

H
CO
I
\
L

H
(D
U

....
K)
I

:
,]

Volume (A3)

H
O

f/

1 iT\1(theory)
TM(eth )

H
(I

 

 

 

14 1 L l l l
0 500 1000 1500 2000 2500 3000

T (K)

 

 

18.0

17.0

Volume (Aa )

...
9‘
O

 

 

 

L

 

14.0 ‘ L
0.0 0.2 0.4 0.6

X

0.8 1.0

 

Figure 4.8: (a) Equilibrium volume V of a Pd atom in bulk Pd as a function of
temperature T. The dashed line is a guide to the eye connecting the calculated data
points. Observed linear expansion of Ag[Pea 58], the neighboring element of Pd, is
shown by the dotted line. (b) Equilibrium volume V per Pd atom in bulk PdH,,
as a function of the hydrogen concentration 2:. The dotted line corresponds to the
observed[Pei 78] expansion coefficient of Pd as a function of 0:. Results of the present
molecular dynamics simulations are given by s in (a) and (b). (From Ref. [Zho 92].
©American Physical Society)

0015345070

' ' = 3001(

 

 

 

 

 

 

 

 

r (A)

°oiééiééiu

92

 

 

 

 

 

 

r (A)

0

T = 205le'

 

 

012345878

r (A)

Figure 4.9: Atomic trajectories in bulk Pd at the temperatures (a) T = 300 K and
(b) T = 1800 K below the melting temperature TM, and (c) at T = 2050 K above
TM. Corresponding pair correlation functions g(r) at T = 300 K, T = 1800 K
and T = 2100 K are shown in (d), (e), and (f), respectively. (From Ref. [Zho 92].
©American Physical Society)

93

an overheated solid or an undercooled liquid, depending on whether the system is
being heated or cooled. The small difference of S 10% between the calculated melting
temperature and the experimental value[Kit 86] of TM = 1827 K is impressively small
in view of the fact that the interaction potentials are based on T = 0 static properties
of Pd. A small positive difference between the calculated and the observed melting
point is also expected due to the finite size of the unit cell and absence of large defects

or a surface.

In the temperature range T < 1500 K, the thermal expansion of the lattice
is nearly linear, corresponding to a thermal (linear) expansion coefficient of a; =
1.7 x 10"5 K". While I could not find any corresponding experimental results for
Pd, this value lies very close to a; = 1.89 x 10‘5 K‘1 which has been observed in
Ag[Pea 58], a neighbor of Pd in the periodic system. At higher temperatures, the
volume-temperature relation shows strong deviations from linearity. This relation-
ship bears information about the Pd-Pd interaction potentials at large interatomic
separations. Above the melting point, the expansion coefficient of the system experi-
ences an abrupt increase above the solid phase value, correlated with a sharp increase

of the interparticle distance at TM.

In order to illustrate the effect of temperature on the dynamical behavior of the
system, I showed the trajectories of individual Pd atoms [projected onto the (100)
plane] in Figs. 4.9(a)-(c) at different temperatures. At T = 300 K, all atoms appear
to be pinned to their equilibrium site and show negligible fluctuations about this po-
sition, reflected in a very small Debye-Waller factor. At T = 1800 K, the fluctuations
of Pd atoms about their equilibrium sites are already quite appreciable. The size of
these fluctuations is a considerable fraction of the lattice constant, which, according

to Lindemann’s criterion, is indicative of the proximity to the melting point. A statis-

94

tically negligible fraction of atoms is also seen to diffuse faraway from their site, but
the crystalline order still exists. Above the melting point, all atoms diffuse relatively
freely throughout the sample, as shown in Fig. 4.9(c), and the long range order is

destroyed.

A more quantitative measure of the crystalline order is the pair correlation func-
tion g(r) which is shown in Figs. 4.9(d)—(f) for the above temperature values. Clearly,
the validity of conclusions related to the long-range order in the crystal is limited by
the finite size of the simulation box. At T = 300 K, g(r) consists of a set of sharp
peaks which will result in a sharp diffraction pattern. At T = 1800 K, g(r) still shows
a substantial amount of structure even for large interatomic distances r, indicative
of long-range order, but the peaks are smeared out and begin to overlap. Above the
melting transition, the most characteristic feature in the pair correlation function is a
peak at the nearest-neighbor distance. For larger interatomic distances, g(r) contains
almost no information about the atomic structure and approaches the constant value
g(r) = 1 very rapidly. This would justify a treatment of the liquid as a continuous

medium beyond the nearest-neighbor distance.

As I will discuss later on in connection with fracture, interesting details about the
cohesion of the system under different conditions can be learned from the distribution
of atomic binding energies. In Figs. 4.10(a) and (b), I compare the binding energy
distribution for bulk Pd atoms at T = 300 K and T = 2100 K, just above the
melting point. In absence of applied tensile stress, my simulations indicate that
binding energies in the system can be characterized by a rather featureless single-
peaked distribution function indicating that most atoms have an indistinguishable
environment. The small width of this peak at T = 300 K, shown in Fig. 4.10(a),

reflects an almost perfect crystalline order at this temperature. Above the melting

95

 

 

 

 

 

 

 

 

 

 

 

T
? 25 :0
‘5' 20(3) ‘80 f-3
'E 33.
€ 15- .50 E
‘3 10 ' -40 ‘5.
3: :-
‘3 5' -2o .:
s ] E
'3', 0 1 0 1
:. 20 25 3.0 3.5 40 31/
Eeoh(ev)
1
3 10 ,- 100
_., 1b) -
.3 81' ‘1‘ .80 E
d ' . ”Q
2 I
E 6- «60 E-
.3. : ' é :
1 ‘r '1' 1‘0 5
3: 'o 1 5'
= .‘ :_-_
:5 2 . I. 7'20 E—
3, 0, #0 J:
‘5. 2.0 . . 4.0 3
T 0
:1 (C) 100
g 8) .30 E-
‘3 a:
:9. s- ‘00 3'-
:53 2
‘3 4* «10 "P.
.. s-
: w
3 2» +20 9,:
0 or a '3
'5. 2.0 2.5 3.0 3.5 4.0 3

Each (0")

Figure 4.10: Distribution of atomic binding energies Eco). in bulk Pd at (a) T = 300 K
and (b) T = 2100 K at zero pressure. (c) Corresponding data for Pd at T = 300 K
at the point of critical uniaxial tensile stress pc for fracture, discussed in Section 4.3.2
The solid lines give the probability distribution, the dashed lines show the integrated
probability. (From Ref. [Zho 92]. ©American Physical Society)

96

point, this peak is strongly broadened and shifted towards smaller binding energies,
as shown in Fig. 4.10(b). The former eflect comes from the large variety of binding
sites in the liquid, the latter one reflects the loss of cohesion mainly due to a uniform

lattice expansion.

4.3.2 Mechanical stability of H-free Pd under tensile stress

One of the most challenging problems related to the mechanical stability of metals is to
obtain a quantitative understanding of the fracture process due to uniaxially applied
tensile stress. I address this problem by studying the deformation of a metal block
(the simulation box) under a uniaxial load, as shown schematically in Fig. 4.11(a). I
expect the length 2 of the metal block first to increase monotonically with increasing
load. Once a critical value of the tensile stress pc is reached, the material can no longer
support the load and breaks into parts. In the following, I describe the molecular
dynamics simulation of the elastic and plastic deformations in the Pd metal block
which is exposed to increasing uniaxial tensile stress, as a function of temperature

and - in the following subsection - as a function of hydrogen concentration.

In the present MD simulation, I consider a canonical (TtN) ensemble of 500 Pd
atoms in a simulation box which has cubic shape at zero applied stress, but which
can vary its shape freely in the case of anisotropic pressure. For a given temperature
of the heat bath, I determine the shape changes [especially the elongation of the cell
Az, see Fig. 4.11(a)] in response to uniaxial tensile stress 'which is described by the

stress matrix

(4.6)

GOO
”606

0
2= 0
0

Results for the elongation A2 of the MD unit cell as a function of p are shown in

 

 

 

97

 

 

 

 

 

 

 

 

 

 

.0 A A A A A
0.0 2.5 5.0 7.5 10.0 12.5150

 

p(GPa)
' - 5". r3200
(0)1151 - rte-01‘
a'osrrz =0. 25
10.0’ t v ‘

 

 

 

 

0. A . L A .
00.0 2.8 8.0 7.8 10.0 12.8 18.0
9(GPn)

Figure 4.11: (a) Schematic picture of Pd deformations under tensile stress, showing
the length 2 and deformation A2 of the unit cell. (b) Deformation A2 of bulk Pd as
a function of the external tensile stress p, for different temperatures. (c) Deformation
A2 of hydrogen loaded bulk Pd at T = 300 K, as a function of the external tensile
stress p, for different hydrogen concentrations. The molecular dynamics results are
given by the data points. The lines are guides to the eye. (From Ref. [Zho 92].
©American Physical Society)

 

 

98

Fig. 4.11(b) for temperatures between 77 K and 1800 K. In this simulation, I increased
the uniaxial stress p in finite steps of initially Ap = 2.0 GPa from zero to a load just
below the point of fracture. From there on, I decreased the steps Ap to 0.5 — 1 GPa
in order to increase the accuracy of the calculated critical tensile stress pc. At each
value of p, the system has been allowed typically 5,000 time steps (or 10-113) to
equilibrate. I observed that the equilibration took longer near the point of fracture
and extended the simulation accordingly. The equilibrium shape of the MD unit cell

has been obtained by averaging over the last 5,000 time steps.

The first important result of this simulation is the order of magnitude for the criti-
cal tensile stress pc to initiate fracture, typically a few GPa. As shown in Fig. 4.11(b),
pc decreases with increasing temperature, from 11 GPa at T = 77 K to 1 GPa at
T = 1800 K. On the other hand, I find an increase in the plasticity, given by Bz/ap,
with increasing temperature at constant load. I conclude that the Young’s modulus Y,
defined by Y = 6p/(92, decreases as the temperature rises. Both effects indicate that
the material becomes softer and easily deformable with increasing temperature, which
agrees with the everyday experience. Microscopically, this softening corresponds to
the increased probability of activated atomic diffusion leading to a new equilibrium
geometry (plastic deformations to a “thin wire” and fracture under excessive uniaxial
load). In the elastic region p << pc, the Az-p relationship is nearly linear for all tem-
peratures. This Hooke’s law type behavior is expected based on the MBA interaction

potential which is dominated by harmonic terms close to the equilibrium.

It is interesting to note that the above MD simulations, performed under uniaxial
stress, contain the information about the bulk modulus B which describes the elastic

response to isotropic pressure. For a cubic crystal, the elastic response to a very small

 

99

applied uniaxial stress p is

p = 01161 + 201262 (4.7)

along the direction of the load. Here, 61,62 are the strain components along the
load direction and perpendicular to the load direction, respectively, and C11, C12 are

elastic stiffness constants. No stress occurs perpendicular to the load direction [see

Fig. 4.11(a)], so

0 = 0126] + 01162 + 01262 . (4.8)

From Eqs. (4.7) and (4.8) one can derive

V(p) - W? = 0)

 

p = (011+ 2012) V(p = 0) (4.9)
Using B = (1/3)(C11 + 2C12) for a cubic crystal leads to
-1
B = 1V (a!) . (4.10)
3 6p F0

As shown in Fig. 4.12(a), the present results for the absolute value of B and the
temperature dependence of B are in good agreement with the experimental data
of Ref. [Lan 79]. At low temperatures, the bulk modulus is found to be essentially
independent of temperature. With rising temperature, however, the present results

indicate a strong decrease of B.

I found it instructive to inspect the distribution of binding energies for a signature
of atomic fracture at very large tensile loads. In Fig. 4.10(c), I show the distribution
at the point of critical tensile stress pc, for a developed fracture. As compared to
the stress-free situation shown in Figs. 4.10(a) and (b), the distribution of binding

energies shows several distinctive peaks. The lowest binding energies correspond to

100

 

 

 

 

200 p. a
(a) ...: ----------------- .
160 ' ~ ”;~ ~:
140 ~ xx] 7
A 120 " “\ 7
(U ‘s l
0- 100 ~ ‘~ 4
8 \ I
m 80 " “x ‘i
60 " “9, '1
40 r i‘:

20 r

o L L 1
0 500 1000 1500 2000
T00

 

200
(b) 190
180
170
180
180
140
130
120
110 '

loo 1 1m L 1m 1 l L l l
0.0 0.2 0.4 0.8 0.8 1.0

 

 

 

 

Figure 4.12: (a) Temperature dependence of the bulk modulus B of Pd. The dashed
line is a guide to the eye connecting the calculated points, given by e. The dotted line
shows the experimental data of Ref. [Lan 79]. (b) Dependence of the bulk modulus of
PdH, on the hydrogen concentration 2:. (From Ref. [Zho 92]. ©American Physical

Society)

 

101

sites at the surface of the crack. The highest binding energies, same as in the stress-

free sample shown in Fig. 4.10(a), correspond to bulk sites in the intact fragments.

4.3.3 Mechanical stability of H loaded Pd under tensile
stress

Hydrogen is well known to dissociate at transition metal surfaces and to penetrate
easily into the bulk metal, releasing the heat of hydride formation in many systems
such as Pd[Ale 78]. This process makes such metals an ideal medium for hydrogen
storage. On the other hand, the presence of hydrogen is known to have an adverse
effect on mechanical properties of metals, specifically facilitating the formation of
cracks under tensile stress[Bir 79]. In order to obtain a microscopic understanding
of the processes associated with hydrogen-assisted crack formation, I simulated the
response of hydrogen loaded Pd to uniaxial tensile stress in a molecular dynamics

calculation.

I have chosen a cubic simulation box containing 500 Pd atoms as the initial MD
unit cell, and occupied the 500 octahedral interstitial sites in the lattice at random
with hydrogen atoms. I have considered three different H concentrations in Pde,
namely :1: = 0.1, :1: = 0.25, and :c = 0.6, and performed all simulations at room
temperature, T = 300 K. At each H concentration, I first let the system equilibrate

over a period of more than 30,000 time steps (corresponding to 1.5 x 10‘113).

First, I study the volume changes due to hydrogen at zero pressure. The free vol-

ume V(z) of Pd atoms in PdH,” shown in Fig. 4.8(b), has been determined for each

H concentration a: by averaging over 10,000 time steps after reaching the equilibrium.

I found the free volume to increase almost linearly with increasing hydrogen concen-

tration x, which agrees with the observed relation[Pei 78] V(x) = V(0)(1 + 0.192;).

 

 

102

In order to understand the effect of dissolved hydrogen on the mechanical stability
of Pd, I studied the response of the system to uniaxial tensile stress using molecular
dynamics. As in the hydrogen-free samples, the system was first allowed to equilibrate
under fixed external tensile stress p and constant temperature of the heat bath T =
300 K. Then, the stress-induced elongation A2 of the simulation box was determined
by averaging over 10,000 time steps (corresponding to 5 x 10’125). The results,
presented in Fig. 4.11(c), show an almost linear relationship between p and A2. I
found the slope of the Az(p) curves to be almost independent of a: at low values
of p, indicating that changes of the hydrogen concentration have little effect on the
Young’s modulus Y. These results also indicate that the critical tensile stress pc
for the onset of fracture, corresponding to the “end points” of the Az(p) curves,
decreases with increasing hydrogen concentration. I found that hydrogen can reduce
the critical tensile stress for fracture pc substantially when compared to the hydrogen-
free system, but that the order of magnitude of pc in the different systems is the
same. This hydrogen-induced reduction of the mechanical strength is sometimes
called “hydrogen embrittlement”[Bir 79]. As I will discuss later on, the microscopic
results indicate that this is a misnomer; I find hydrogen to enhance the ductility and

plasticity of the metal matrix locally, thereby weakening the structure as a whole.

One can also use the MD results for uniaxial tensile stress to estimate the bulk
modulus of the system, following the procedure outlined in the previous subsection.
The results, presented in Fig. 4.12(b), indicate that the bulk modulus B decreases
strongly in the presence of hydrogen. Another quantity of interest is the Poisson’s
ratio )1 which is the ratio of the unit cell deformations along the direction of the

applied load, A2, and perpendicular to it, Am. In a cubic system, )1 relates the

 

103
Young’s and the bulk modulus as B = Y/ (3 — 6p), which defines p as

L
63'

11 = (4.11)

NIH

My above results for Y(2:) and B(:r) indicate that the Poisson’s ratio decreases

strongly with increasing hydrogen concentration in the metal.

As in the hydrogen—free system, I investigated the distribution of binding energies
for a signature of atomic fracture at different tensile loads. My results for the PdHo,25
system are summarized in Fig. 4.13(a) for zero tensile stress and in Fig. 4.13(b) for
critical tensile stress. The structure of the binding energy distribution in the stress-
free case reflects the distribution of inequivalent Pd atoms with 0—6 hydrogen nearest
neighbors. The relatively featureless distribution of binding energies at the point of
critical tensile stress p = pc, displayed in Fig. 4.13(b), is in sharp contrast to the
hydrogen-free case shown in Fig. 4.10(c), and is more reminiscent of the results for
molten Pd, shown in Fig. 4.10(b). The absence of distinct fracture-related features in
the binding energy distribution in Fig. 4.13(b) indicates that all Pd atoms reside in a
relatively homogeneous atomic environment. This environment is close to the molten
system; it has an amorphous structure and can easily be plastically deformed. I
conclude that increased hydrogen concentration has a similar effect on the mechanical

properties of Pd as a temperature increase, namely enhanced ductility and plasticity.

An independent microscopic signature of fracture on the atomic scale can be of-
ten found in the distribution of Wigner-Seitz volumes associated with the individual
atoms in the crystal. Statistical presence of large atomic volumes indicates the oc-
currence of fracture with no additional assumptions about the number, position or
morphology of one or more simultaneous cracks. This information is displayed in

Fig. 4.14 for Pd and H atoms in PdHogs at p = pc. The volumes of Pd atoms, given

 

104

 

 

 

 

 

 

 

 

   

 

 

 

 

'> 25]
:7 l (8.) :°
‘5‘ 20] z-
: [ a-s.
: '3
73 15: 5;:
:2 l f:
“'3 10 . .3
3» 5-
:: l. a:
'3 5 E
_g :2
a 0 ‘<
a. 2.0 2.5 3.0

E:eoh (8V)
T
E 10 (b) I, 100
5 8 . 5 80 §
‘3 : s
g 6- 5 ‘60 g»
72‘ 5 °‘
1 4 L ‘40 E
3: E—
: 2 _ ‘ 0a
3 , 20 ___c-_:
s :-
g o - e . 0 7‘
ca. 2.0 2.5 3.0 3.5 4.0

Ecoh (8V)

Figure 4.13: Distribution of binding energies of Pd atoms in bulk PdHo_25 at T =
300 K (a) at zero pressure and (b) at the point of critical uniaxial tensile stress
pc [compare with Fig. 4.11(c)]. The solid lines give the probability distribution,
the dashed lines show the integrated probability. (From Ref. [Zho 92]. ©American

Physical Society)

 

 

 

105

.0
on
0
fl
1

 

.0
10
CI

0.15

0.10

0.05

 

 

 

 

0.00

 

Volume distribution (A’s)

 

Figure 4.14: Distribution of atomic volumes in bulk PdHo,25 at T = 300 K. The solid
line shows the distribution of atomic volumes associated with Pd atoms, obtained
using the Wigner-Seitz cell construction. The dashed line is obtained by first deter-
mining, which Pd sites are associated with each of the H atoms in the lattice, and
displaying the distribution of these Pd volumes. The dotted line gives the ratio of the
values given by the dashed and the solid lines, divided by the volume and multiplied
by 3. (From Ref. [Zho 92]. ©American Physical Society)

 

106

by the solid line, show a distribution which is characterized by a sharp peak near
15 A3, corresponding to atoms in bulk-like environment, and a wide structureless tail
towards larger volumes, associated with atoms near the crack surface. This finding
again confirms the previous conclusion that Pd atoms have no preferential binding
arrangement in the hydrogen loaded sample, which shows a plastic behavior under

critical tensile stress.

The preferential hydrogen sites in this structure are determined in the following
way. I enlarged the original Wigner-Seitz volumes of Pd atoms by a factor of 1.5 in
each direction and for each hydrogen atom in the crystal, I generated a list of Pd sites
which contain this atom] in their enlarged unit cell. This definition allows for more
than one Pd site to be associated with a given hydrogen atom. Next, I combined
the lists of Pd sites associated with each H atom and plotted the distribution of the
corresponding Pd Wigner-Seitz volumes. The results are given by the dashed line
in Fig. 4.14. A comparison with the solid line for the Pd atoms shows no strong
preference of hydrogen atoms for specific sites in the metal structure. The dotted line
in Fig. 4.14 represents the probability that a Pd atom, characterized by its atomic
volume, is likely to have one or more hydrogen atoms as its closest neighbors. In
case that there would be no preferential sites for hydrogen atoms, this curve should
be flat. These results indicate a strong preference of hydrogen atoms for sites near
highly coordinated Pd atoms in the bulk. From the position of the peak in the dotted
curve, which lies at a slightly larger volume than that of bulk Pd atoms, one can infer
that hydrogen atoms are more likely to occupy subsurface sites or sites close to the

crack tip than sites in the bulk of Pd.

107

4.3.4 Discussion

The above molecular dynamics calculations, based on the N osé and Rahman-Parrinello
formalism, suggest a useful and consistent picture of the atomic-scale processes which
occur in PdH, at different temperatures and hydrogen concentrations as, specifically
in response to large uniaxial tensile stress. Even though the interaction potentials
are based on static ab initio calculations at T = 0, the finite temperature results of
these simulations are in good agreement with experimental data. This increases my
confidence that the MBA Hamiltonian describes the interactions in the Pd—H system
correctly to a large degree. Since this Hamiltonian has a solid theoretical background,
it can provide microscopic insight into the nature of interatomic interaction under dif-

ferent conditions.

These simulations show that the introduction of hydrogen into bulk Pd at room
temperature has a similar effect on the structural properties and elastic behavior
as a temperature increase in the hydrogen-free metal. Increased hydrogen loading
and increased temperature both increase the free volume linearly, decrease the bulk
modulus, and reduce the critical tensile stress for fracture. It is plausible to some
degree that the presence of hydrogen simulates a temperature increase, since the
light H atoms have a much faster dynamics that Pd atoms and can easily excite Pd
vibrations in elastic collisions. This effect is enhanced by the fact that hydrogen
atoms reduce the bonding strength between Pd atoms and soften the vibrational
modes of the Pd lattice[Zho 91a]. I found that hydrogen-induced changes of structural
properties and elastic response become more pronounced with increasing hydrogen
concentration, as shown in Figs. 4.8(b), 4.12(b) and 4.11(c). At the point of fracture, I
found that the presence of hydrogen enhances the plasticity of the system significantly,

causing hydrogen-assisted “melting” even at T = 300 K, which is associated with a

108

substantial diffusion of Pd atoms. I found it useful to‘compare the present results for
hydrogen loaded Pd under critical tensile stress to hydrogen-free Pd at the melting
point in an animated video movie, and found Pd diffusion at the crack surface of the
hydrogen loaded system to be comparable with the atomic diffusion in melting Pd

metal.

As mentioned above, a fundamental difference between the effect of hydrogen
loading and temperature increase lies in the fact that hydrogen modifies the interac-
tion between Pd atoms. This difference is most obvious in the response to uniaxial
tensile stress, shown in Figs. 4.11(b) and 4.11(c). In the hydrogen-free system, one
can observe both the bulk and the Young’s modulus to decrease with increasing tem-
perature, while the Poisson’s ratio p is nearly constant, consistent with Eq. (4.11).
On the other hand, increasing the hydrogen concentration in Pd at a constant tem-
perature T = 300 K still causes the bulk modulus to decrease, but has little effect on
the Young’s modulus. In this case, the Poisson’s ratio decreases with increasing hy-
drogen loading, indicating an increasing resistance against shape deformations. This
effect could also assist in the initial formation of cracks. Based on the results shown in
Fig. 4.11(c), I find that the critical value of tensile stress pc drops at an increasing rate
with increasing hydrogen concentration, indicating a continuous ductility increase in

the system with increasing number of H atoms in the vicinity of the Pd metal bonds.

The microscopic origin of the ductility increase of Pd in presence of hydrogen is
the rehybridization of Pd orbitals in the hydrogen loaded metal. As was found in
a previous ab initio calculation of atomic binding in bulk Pd and PdH[Sun 89], the
hydrogen-induced binding changes of Pd stem mainly from a filling of antibonding
states in the Pd4d band, accompanied by a depletion of the partly filled Pd5s band

and a small charge transfer towards hydrogen. These effects are considered, albeit

...“,

 

109

in a very approximate way, in the many-body alloy Hamiltonian. Unlike in models
based on pairwise interactions, this Hamiltonian does address the different electronic
hopping processes to neighboring atoms from the point of view of electronic band
formation in the attractive part of the total energy. Hence the binding energy of a Pd
atom is not simply proportional to the number of nearest neighbors, but depends in
a more complex way on the hybridization with the neighboring atoms and the band
filling, which appears to be essential for the understanding of bonding changes in this

system.

The calculated values for the critical tensile stress are about one order of mag-
nitude too high when compared to experimental data for single- and polycrystalline
samples. There are two reasons for this overestimate of pc. First, it is impossible
to study the dynamics of the fracture process in a molecular dynamics calculation
for realistic time scales and large systems. The time spans presently accessible by
such simulations fall at least ten orders of magnitude short of a realistic time for the
formation of a crack, which is typically seconds. Second, in a realistic system, the
fracture process is assisted and proceeds by dislocation motion in the crystal. The
presently used unit cell containing 500 Pd atoms, which is large for MD standards, is

clearly too small to show the spontaneous formation and motion of dislocations.

I also tried to address the effect of dislocations in a somewhat artificial way,
namely removing a single atom from the unit cell, thus creating a defect site and a
possible seed for a dislocation or a crack. The corresponding MD simulation did not
show a reduction of pc when compared to the initially “perfect” systems. This is not
surprising since the typical atomic fluctuations near the point of fracture are large

and comparable to the size of a single atomic vacancy.

In hydrogen loaded Pd, the local decohesion and structural relaxations are medi-

 

 

110

ated by hydrogen which has a very large diffusion constant. I studied the interplay
between the time scales for structural changes and the diffusion of hydrogen by con—
sidering isotope substitution. I found that replacing H by D atoms in PdH,c has no
effect on pc, probably due to the large difference between the time scales for hydrogen
motion and structural relaxation, and possibly also the incoherent motion of hydrogen

atoms in the metal which averages out local changes of elastic properties.

The above microscopic results for the elastic response of hydrogen-loaded Pd
to uniaxial tensile stress strongly support one of the previously postulated mecha-
nisms for “hydrogen embrittlement”, namely the Hydrogen Enhanced Local Plasticity
(HELP) mechanism[Bir ]. This mechanism, which has been used to interpret exper—
imental data, postulates that hydrogen concentrates preferentially near the tip of a
starting crack. Hydrogen subsequently locally softens the metal matrix in the vicinity
of the crack tip, which leads to an increase in the velocity of dislocation motion. This
process can lead to a softening over microns over very short time scales. The system
will become microscopically ductile, but will appear as brittle on macroscopic length

scales.

While the present calculations have been performed for a specific system, namely
hydrogen loaded Pd, I expect that the effect of hydrogen on the stability of other fcc
metals will be qualitatively the same. The situation in bcc metals (such as Fe) may
be somewhat different, since hydrogen is observed to stiffen rather than soften these

structures[Bir ]. These questions are presently being addressed in corresponding MD

simulations[Zho 93b].

 

 

 

111

4.4 Summary and conclusions

In conclusion, I applied the Many-Body Alloy (MBA) Hamiltonian to Pd—H systems.
All parameters have been obtained from ab initio Density Functional calculations,
with no adjustable parameters for surface properties. I tested this Hamiltonian first
and calculated the equilibrium structure and binding energy of bulk Pd and PdH,
as well as H-free and H covered Pd(001) and Pd(110) surfaces. I found the calcu-
lated results to be in good agreement with experimental data and results of ab initio

calculations where available.

Next, I studied the effect of hydrogen on the vibration spectra of bulk Pd and
the Pd(001) and Pd(110) surfaces, by constructing the dynamical matrix based on
the MBA Hamiltonian. I found that in the bulk systems, hydrogen softens the Pd
vibration modes, as seen in the comparison of bulk Pd and PdH phonon spectra. The
results for the clean and H covered (001) and (110) surfaces of Pd are not as clear-cut.
I found the most pronounced eflect of hydrogen coverage to be the softening of the
surface Rayleigh mode with out-of-plane vibration amplitudes on the topmost layer.
Other surface modes, such as in-plane vibrations of the topmost layer, are affected to

a lesser degree or occur at higher vibration frequencies in the presence of hydrogen.

Finally, I have used the N osé and Rahman-Parrinello molecular dynamics formal-
ism to study the equilibrium structure and elastic properties of bulk Pd as a function
of temperature and hydrogen concentration. I have used this formalism first to predict
the elastic constants, thermal expansion and melting temperature of hydrogen-free
and hydrogen loaded bulk Pd. I found my results to be in good agreement with

available experimental data.

Introducing uniaxial tensile stress as an independent variable into this formalism

 

 

112

has enabled me also to study the elastic deformations as a function of the applied
load at different temperatures and hydrogen concentrations. At small applied loads,
I found that the bulk and the Young’s moduli decrease with increasing temperature
in hydrogen-free bulk Pd. Increased hydrogen concentration at constant temperature
has a very similar effect as the temperature increase in the hydrogen-free metal: The
system gets softer, which is reflected in a decreased bulk modulus. While hydrogen
softens the Pd-Pd bonds, its presence does not affect the Young’s modulus. Conse-
quently, the Poisson’s ratio decreases with increasing hydrogen loading, indicating an
increasing resistance towards shape deformations. This behavior might assist in the

formation of cracks.

At large values of the uniaxial tensile stress, one can observe the onset of crack for-
mation. I find that the critical tensile stress for fracture decreases both with increasing
temperature and increasing hydrogen concentration. Near the point of fracture, how-
ever, the elastic response of hydrogen-free and hydrogen-loaded Pd is vastly different.
Following the fracture, Pd atoms can be found in well-defined “crystalline” sites in
the fragments or at the crack surface for the Pd system which is free of hydrogen.
In hydrogen-loaded Pd, the metal structure can be called amorphous at the point
of fracture. One can find a broad distribution of Pd sites in this system which can
easily be deformed plastically. I conclude that the hydrogen-induced reduction of the
mechanical stability of Pd (and likely also other metals) originates from an increased
ductility and plasticity in parts of the sample with a large hydrogen concentration,
such as regions near grain boundaries and dislocations. These conclusions agree with
one of the previously postulated mechanisms for “hydrogen embrittlement”, namely
the Hydrogen Enhanced Local Plasticity (HELP) mechanism[Bir ]. More detailed

studies will be necessary to confirm this behavior also in bcc metals (such as iron). A

 

113

comparison with corresponding experimental data, stemming from atomic resolution

studies of single crystals would be highly desirable.

 

Bibliography

[Abr 89] F.F. Abraham and LP. Batra, Surf. Sci. 209, L125 (1989).
[Alb 87] T.R. Albrecht and CF. Quate, J. Appl. Phys. 62, 2599 (1987).

[Ale 78] G. Alefeld and J. Valkl, Eds., Hydrogen in Metals I and II, Vols. 28 and 29

of Topics in Applied Physics (Springer—Verlag, Berlin, 1978).

[A1190] M.P. Allen and DJ. Tildesley, Computer Simulation of Liquids (Oxford
Press, New York, 1990).

[And ] M. Anders and C. Heiden (submitted for publication).
[And 80] H.C. Andersen, J. Chem. Phys. 72, 2384 (1980).

[Arn 66] RD. Arnell, J.W. Midgley, and D.G. Teer, Proc. Inst. Mech. Eng. 179, 115
(1966).

[Bac 82] GB. Bachelet, D.R. Hamann, and M. Shliiter, Phys. Rev. B 26, 4199 (1982).

[Bar 85] C.J.Barnes, M.Q. Ding, M. Lindroos, R.D. Diehl, and DA. King, Surf. Sci.
162, 59 (1985).

[Bas 89] R. Bastasz, T.E. Felter, and WP. Ellis, Phys. Rev. Lett. 63, 558 (1989).

[Bat 88] LP. Batra and S. Ciraci, J. Vac. Sci. Technol. A 6, 313 (1988).

114

 

115

[Bea 72] CD. Beachem, Metall. Trans. 3, 437 (1972).
[Beh 80] R.J. Behm, K. Christmann and G. Ertl, Surf. Sci. 99, 320 (1980).
[Bes 87] F. Besenbacher, I. Stensgaard and K. Mortensen, Surf. Sci. 191, 288 (1987).

[Bin 82] G. Binnig, H. Rohrer, Ch. Gerber and E. Weibel, Phys. Rev. Lett. 49, 57
(1982).

[Bin 86] G. Binnig, C.F. Quate and Ch. Gerber, Phys. Rev. Lett. 56, 930 (1986),
and Appl. Phys. Lett. 40, 178 (1982).

[Bin 87] G. Binnig, C. Gerber, E. Stoll, T. R. Albrecht, C. F. Quate, Europhys. lett.
3, 1281 (1987).

[Bir 73] H. K. Birnbaum, M. Grossbeck, and S. Gahr, in: Hydrogen in Metals, M.
Bernstein and A. Thompson, eds., A.S.M. Metals Park, Ohio, 1973, p. 303.

[Bir 79] H.K. Birnbaum, in: Environmentally Sensitive Fracture of Engineering Ma-

terials, Z. A. Foroulis, eds., T.M.S. New York, 1979, p. 326.
[Bir 84] H.K. Birnbaum, J. Less Common Metals 103, 31 (1984).

[Bir ] Howard K. Birnbaum, private communication.

[Bob 90] K.P. Bohnen, Th. Rodach, and KM. H0, in The Structure of Surfaces [11,
edited by M.A. Van Hove, K. Takanayagi, and X. Xie (Springer-Verlag, Hei-

delberg, 1990), and references cited therein.

[Bor 84] V. Bortolani, G. Santoro, U. Harten, and JP. Tonnies, Surf. Sci. 148, 82
(1984).

 

 

116

[Bor 85] V. Bortolani, A. Franchini, F. Nizzoli, and G. Santoro, Surf. Sci. 152/153,
811 (1985).

[Bre 89] D.W. Brenner, Phys. Rev. Lett. 63, 1022 (1989).
[Cal 72] Langlinais and J. Callaway, Phys. Rev. B 5 124 (1972).
[Cal 84] J. Callaway and NH. March, Solid State Phys. 38, 135 (1984).

[Car 90] A.E. Carlsson, in Solid State Physics, Vol. 43, Academic Press (New York,
1990), p. 1.

[Cat 83] M.G. Cattania, K. Christmann, V. Penka and G. Ertl, Gazz. Chim. Ital.
113,433 (1983).

[Cep 80] D.M. Ceperley and B.J. Alder, Phys. Rev. Lett. 45 566 (1980).
[Cha 73] DJ. Chadi and M.L. Cohen, Phys. Rev. B 8, 5747 (1973).
[Cha 86a] C.T. Chan, D. Vanderbilt, and S. G. Louie, Phys. Rev. B 33, 2455 (1986).

[Cha 86b] C.T. Chan, D. Vanderbilt, S.G. Louie and JR. Chelikowsky, Phys. Rev.
B33, 7941 (1986).

[Che 84] J. R. Chelikowsky and S. G. Louie, Phys. Rev. B 29, 3470 (1984).

 

[Cir 90] S. Ciraci, A. Baratoff and LP. Batra, Phys. Rev. B 41, 2763 (1990).
[Dav 83] H.L. Davis, J.R. Noonan, Surf. Sci. 126, 245 (1983).

[Daw 84] MS. Daw and M.I. Baskes, Phys. Rev. B 29, 6443 (1984).

[Daw 89] MS. Daw, Phys. Rev. B39, 7441 (1989).

[Doa 83] RB. Doak, U. Harten and JP. ’Toennis, Phys. Rev. Lett. 51, 578 (1983).

 

117

[Dre 90] RM. Dreizler and E.K.U. Gross, Density Functional Theory (Springer-
Verlag, 1990).

[Fah 86] S. Fahy, S.G. Louie and M.L. Cohen, Phys. Rev. B 34, 1191 (1986).
[Foi 85] SM. Foiles, Phys. Rev. B32, 3409 (1985).

[Foi 86] SM. Foiles, M.I. Baskes, and MS. Daw, Phys. Rev. B 33, 7983 (1986).
[Foi 87] S.M. Foiles, Surf. Sci. 191, L779 (1987).

[Foi 89] SM. Foiles and JB. Adams, Phys. Rev. B40, 5909 (1989).

[Gea ] C.W. Gear, private communication.

[Con 89] S. Gould, K. Burke, P.K. Hansma, Phys. Rev. B 40 5363 (1989).

[Gru 89] M. Grunze and H.J. Kreuzer, Eds., Adhesion and Friction, Springer Series

in Surface Sciences, Vol. 17 (Springer-Verlag, Berlin 1989).

[Hal 88] BM. Hall, D.L. Mills, M.H. Mohamed, and LL. Kesmodel, Phys. Rev. B38,
5856 (1988).

[Ham 79] DR. Hamann, M. Schliiter and C.Chiang, Phys. Rev. Lett. 43, 1494 (1979).
[Har 82] B. Harmon, W. Weber, and D. R. Hamann. Phys. Rev. B 25, 1109 (1982).

[Har 85] U. Harten, J .P. Toennis, Ch. W611 and G. Zhang, Phys. Rev. Lett. 55, 2308
(1985).

[He 88] J.—W. He, D.A. Harrington, K. Griffiths, and RR. Norton, Surf. Sci. 198,
413 (1988).

[Hed 71] L. Hedin and B.J. Lundqvist, J. Phys. C4, 2064 (1971).

 

118

[Ho 86] K.-M. Ho and K.P. Bohnen, Phys. Rev. Lett. 56, 934 (1986).
[Hoh 64] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).

[Her 92] C. Horie and H. Miyazaki (submitted for publication).

[Iba 87] H. Ibach, J. Vaccum Sci. Technol. A5 419 (1987).

[Joh 72] RA. Johnson, Phys. Rev. B6, 2094 (1972).

[Jon 89] RD. Jones and O. Gunnarsson, Rev. Mod. Phys. 61, 689 (1989).

[Kit 86] Ch. Kittel, Introduction to Solid State Physics, 6th edition (John Wiley, New
York, 1986).

[Koh 65] W. Kohn and L.J. Sham, Phys. Rev. 140 A1133 (1965).
[Kus 90] D. Kusnezov, A. Bulgac and W. Bauer, Ann. Phys. 204, 155 (1990).

[Laf 71] Lafon, R.C. Chaney, and C.C. Lin, in Computational Methods in Band The-
ory, edt. P.M. Marcus, D.F. Janak, and AR. Williams (Plenum, New York,
1971), p. 284.

[Lah 87] A.M. Lahee, J.P. Toennies and Ch. W611, Surf. Sci. 191, 529 (1987).

[Lan 60] L.D. Landau and EM. Lifshitz, Course of Theoretical Physics, Volume 1:
Mechanics (Pergamon, Oxford 1960), p. 122.

[Lan 86] L.D. Landau and EM. Lifshitz, Course of Theoretical Physics, Volume 7:
Theory of Elasticity, (Pergamon, Oxford 1986), p. 53.

[Lan 79] Landolt—Bérnstein (new series), edited by K.-—H. Hellwege (Springer—Verlag,
New York, 1979), Vol. III/11, pp. 10, 114.

 

119

[Lan 89a] U. Landman, W.D. Luedtke, A. Nitzan, Surf. Sci. 210, L177 (1989).

[Lan 89b] Uzi Landman, W.D. Luedtke, and M.W. Ribarsky, J. Vac. Sci. Technol. A
7, 2829 (1989).

[Lan 90] U. Landman, W.D. Luedtke, N .A. Burnham, and R.J. Colton, Science 248,
454 (1990).

[Lee 89] S. Lee, H. Miyazaki, S.D. Mahanti and SA. Solin, Phys. Rev. Lett. 62, 3066
(1989). '

[Leh 83] S. Lehwald, J.M. Szeftel, H. Ibach, T.S. Rahman and D.L. Mills, Phys. Rev.
Lett. 50, 518 (1983).

[Leh 87] S. Lehwald, B. Voigtliinder, and H. Ibach, Phys. Rev. B 36, 2446 (1987).
[Lev 79] M. Levy, Proc. Natl. Acad. Sci. (USA) 76, 6062 (1979).

[Lun 83] S. Lundqvist and N .H. March, Theory of the Inhomogeneous Electron Gas
(Plenum, 1983).

'[Luo 88] N. Luo, W. Xu, and S.C. Shen, Solid State Commun. 67, 837 (1988).
[Lyn 86] S.P. Lynch, J. Mat. Sci. 21, 692 (1986).

[Mah 90] G.D. Mahan and KR. Subbaswamy, Local Density Theory of Polarizability
(Plenum, 1990).

[Mam 86] H.J. Mamin, E. Ganz, D.W. Abraham, R.E. Thomson, and J. Clarke,
Phys. Rev. B 34, 9015 (1986).

[Mar 63] A.A. Maradudin, E.W. Montroll, and G.H. Weiss, Theory of Lattice Dy-

namics in the Harmonic Approximation (Academic Press, New York, 1963).

 

120

[Mat 87] C. Mathew Mate, Gary M. McClelland, Ragnar Erlandsson and Shirley
Chiang, Phys. Rev. Lett. 59, 1942 (1987).

[McC 89] G.M. McClelland in Adhension and friction, edited by M. Grunze and
H. j. Kreuzer, Springer Series in Surface Science, Vol. 17 (Springer-Verlag,
Berlin) 1989.

[Mey 88] E. Meyer, H. Heinzelmann, P. Griitter, Th. Jung, Th. Weiskopf, H.-R.
Hidber, R. Lapka, H. Rudin, and H.-J. Giintherodt, J. Microsc. 152, 269
(1988).

[Mey 90] G. Meyer and N. M. Amer, Appl. Phys. Lett. 56, 2100 (1990).

[Mil 71] A.P. Miller and B.N. Brockhouse, Can. J. Phys. 49,704 (1971).

[Nel 77] J.S. Nelson, E.C. Sowa, M.S. Daw, Phys. Rev. Lett. 61, 1977 (1988).
[Nel 89] J.S. Nelson, M.S. Daw, and EC. Sowa, Phys. Rev. B40, 1465 (1989).
[Nos 84] s. Nosé, J. Chem. Phys. 81, 511(1984); Mod. Phys. 52,255 (1984).
[Nyb 83] C. Nyberg and C.C. Tengstiil, Phys. Rev. Lett. 50, 1680 (1983).
[Ori 77] RA. Oriani and RH. Josephic, Acta Metal]. 25, 979 (1977).

[Ove 92] RM. Overney, E. Meyer, J. Frommer, D. Brodbeck, R. Luthi, L. Howald,
H.-J. Guntherodt, M. Fujihira, H. Takano, and Y. Gotoh, Nature 359, 133
(1992)

[Par 80] M. Parrinello and A. Rahman, Phys. Rev. Lett. 45, 1196 (1980).
[Par 81] M. Parrinello and A. Rahman, J. Appl. Phys. 52, 7182 (1981).

[Par 82] M. Parrinello and A. Rahman, J. Chem. Phys. 76, 2662 (1982).

 

121

[Pea 58] W.B. Pearson, A Handbook of Lattice Spacings and Structures of Metals and
Alloys, Pergamon, New York, 1958.

[Pei 78] H. Peisl in Hydrogen in Metals I, Vol. 28 of Topics in Applied Physics, edited
by G. Alefeld and J. Vélkl (Springer-Verlag, Berlin, 1978), p. 69.

[Pet 87] J.B. Pethica, W.C. Oliver, Physica Scripta T 19, 61 (1987).
[Pra 89] LR. Pratt and J. Eckert, Phys. Rev. B 39, 13170 (1989).
[Ray 84] J.R. Ray and A. Rahman, J. Chem. Phys. 80, 4423 (1984).
[Ray 88] J.R. Ray, Comp. Phys. Rep. 8, 109 (1988).

[Rie 83] K.H. Rieder, M. Baumberger and W. Stocker, Phys. Rev. Lett. 51, 1799
(1983).

[Rie 84] K.H. Rieder and W. Stocker, Surf. Sci. 148, 139 (1984).

[Roc 84] S. M. Rocca, S. Lehwald, H. Ibach and T.R. Rahman, Surf. Sci. 138, L123
(1984).

[Row 74] J.M. Rowe, J.J. Rush, H.C. Smith, M. Mostoller, and HE. Flotow, Phys.
Rev. Lett. 33, 1297 (1974).

[Row 86] J.M. Rowe, J.J. Rush, J.E. Schirber and J.M. Mintz, Phys. Rev. Lett. 57,
2955 (1986).

[Ski 71] J. Skinner, N. Gane, and D. Tabor, Nat. Phys. Sci. 232, 195 (1971).

[Sko 87] M. Skottke, R.J. Behm, G. Ertl, V. Penka, and W. Moritz, J. Chem. Phys.
87,6191 (1987).

122

[Sla 93] M. Slameron, R. Overeny, E. Meyer, M.O. Robbins, P.A. Thompson, C.C.
Crest, J.A. Harrison, C.T. White, R.J. Colton, D.W. Brenner, J. Belak, D.B.
Boercker, I.F. Stowers, U. Landman and W.D. Luedtke, MRS Bulletin, May
1993, p. 15.

[Smo 41] R. Smoluchowski, Phys. Rev. 60, 661 (1941).
[Spa 84] D. Spanjaard and M.C. Desjonqueres, Phys. Rev. B30, 4822 (1984).
[Spr 62] J. Spreadborough, Wear 5, 18 (1962).

[Ste 60] EA. Steigerwald, F.W. Schaller, and AR. Troiano, Trans. Metal]. Soc..
A.I.M.E. 218, 832 (1960).

[Sto 89] P. Stoltze, J .K. N ¢rskov and Uzi Landman, Surf. Sci. Lett. 220, L693 (1989).
[Sun 89] Z. Sun and D. Tomanek, Phys. Rev. Lett. 63, 59 (1989).

[Tab 84] T. Tabata and H.K. Birnbaum, Scripta Metal]. 18, 231 ( 1984).

[Tar 86] B. Tardy, J.C. Bertolini, CR. Acad. Sci., Ser. 2, 302 813 (1986).

[Ter 83] J. Tersoff and DR. Hamann, Phys. Rev. Lett. 50, 1998 (1983), and Phys.
Rev. B31, 805 (1985).

[Toe 87] J.P. Toennies, J. Vacuum Sci. Technol. A5, 440 (1987); ibid. A2, 1055
(1984).

[Tom 83] D. Tomének, S. Mukherjee and K.H. Bennemann, Phys. Rev. B28, 665
(1983); ibid. 829, 1076 (1984) (E).

[Tom 85a] D. Tomanek, A.A. Aligia and C.A. Balseiro, Phys. Rev. B32, 5051 (1985).

[Tom 85b] D. Tomanek and K.H. Bennemann, Surf. Sci. 163, 503 (1985).

123

[Tom 86a] D. Tomanek, Phys. Lett. 113 A, 445 (1986).
[Tom 86b] D. Tomének, S.C. Louie, and C.T. Chan, Phys. Rev. Lett. 57, 2594 (1986).
[Tom 86c] D. Tomének and S.C. Louie, Phys. Rev. Lett. 57, 2594 (1986).

[Tom 89] D. Tomének, G. Overney, H. Miyazaki, S.D. Mahanti and H.J. Giintherodt,
Phys. Rev. Lett. 63, 876 (1989) and ibid. 63, 1896(E) (1989).

[Tom 91a] D. Tomének, Z. Sun, and S.C. Louie, Phys. Rev. B43, 4699 (1991).
[Tom 91b] D. Tomanek and W. Zhong, Physical Review B (in press).

[Tom 91c] D. Tomanek, W. Zhong, and H. Thomas, Europhys. Lett. 15, 887 (1991).
[Tom] 29] C.A. Tomlinson, Phil. Mag. S. 7, Vol. 7, 905 (1929).

[Wes 69] D.G. Westlake, Trans. Am. Soc. Metals 62, 1000 (1969).

[Wig 38] ER Wigner, Trans. Faraday Soc. 34 678 (1938).

[Wut 86] M. Wuttig, R. Franchy and H. Ibach, Solid State Commun. 57, 445 (1986).

[Yan 91] Liquiu Yang, Talat S. Rahman and Murray S. Daw, Phys. Rev. B44, 13725
(1991).

[Yos 88] J. Yoshinobu, M. Onchi and M. Nishijima, Phys. Rev. B 38, 1520 (1988).

[Zab 89] H. Zabel in Graphite Intercalation Compounds, Topics in Current Physics,
edited by H. Zabel and S.A. Solin (Springer-Verlag, New York, 1989).

[Zho 90] W. Zhong and D. Tomanek, Phys. Rev. Lett. 64, 3054 (1990).
[Zho 91a] W. Zhong, Y.S. Li and D. Tomanek, Phys. Rev. B 44, 13053 (1991).

[Zho 91b] W. Zhong, G. Overney, and D. Tomanek, Europhys. Lett. 15, 49 (1991).

124

[Zho 92] W. Zhong, Y. Cai, and D. Tomanek, Phys. Rev. B 46, 8099 (1992).
[Zho 93a] W. Zhong, G. Overney, and D. Tomének, Phys. Rev. B 47 (1993).

[Zho 93b] W. Zhong and D. Tomének, in preparation.

"‘9199559“