‘13....» .. . 1

6

w 1 . m ‘IL‘
a§s§2§ v.3; .
:uxwmvm... 1m“
, . . . . .. ,3. Jam
, :
5.3.2

. 4' 0N1»

7:.

. Pd

. 4. ix
. . . s
.1 z: ,

antral. x! .

3c

..
I. E.
4

.id...
\hflvwwi
r. s:

.. .y
, . 1 4mm... . fink.m.t.r.m:n$1u
, I .m’varnfi.
. . . . . g It!“
. , I

rhruu. , ,
‘ u~ffl..:...n.fi4h.v .

< 4.!
- 3 ..,
~ 0.1: . It.
.. I. F
1...: i; .3

uxfiuuw;
t... 7

.Q I.)
. «Ii-It
. 131.5
‘ , 5..
13.5.1. . ‘
.19... 33.5.25
:3: If: x

a.
9?;

,2 :5)

k:

4....

.nru‘

u? r
.

, :52...

.3 z...

3.

‘ in
3.33.... r . ..
5......msn’av “waived
11... .2. 23.5.3.
9?: .5.
I
.3}... .-
...{cxsligi .3298.\ d:
a ofini is! a
,itimz). .41. .3. .. vi? ‘2.
A... 3:. F}...:‘2i.: tyre
.r..1.;1\ie 1.7! Amu‘io
.O‘Ok' I, .9!

1
v .5.»

u.” an
15...:

X.

.

52.13.14
A to .
a 2‘3qu i...
'0 in
. !) nu. Ida.
, . , rg‘
{#hvin. .23.}! ‘
rtkdrc f

i: 27).}? l .
DQ301513}... 510...?!
3: i)! o!
I... Wait:
I a

41.]. .7...)
717.1?63

1!

{1 Av?

i... firh. “kidntfid 1'
{1.0.9231 2!:

”90““:th

.u'mu

:}

“I

>1... any.
\(al

“(I It

 

Z .1
l0.

 

 

 

THcSIS

lHlHill!"IllHIUHHIHIJHIIIIWNllllllllllllll

301046 1600

This is to certify that the
dissertation entitled
MORPHOLOGIES OF

METAL SURFACES, FILMS AND CLUSTERS

presented by

Xinhua Yu

has been accepted towards fulfillment
of the requirements for

Ph . D . degree in Physics

 

 

fiZM/W)

Major professor/

Date /{//flI/q?

MS U is an Affirmative Action/Equal Opportunity Institution 0-12771

 

LIBRARY
Mlchlgan State
Unlverslty

 

 

 

PLACE II RETURN Boxwmwombchockouflunywmd.
TO AVOID FINES atom on orbdorodatoduo.

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU loAnNflmutlvo Action/Equal Opponmlty Intuition
W1

k..__fi ,__— _ — 7 A- ,

 

MORPHOLOGIES OF
METAL SURFACES, FILMS AND CLUSTERS

BY

Xinhua Yu

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

1 994

ABSTRACT

MORPHOLOGIES OF METAL FILMS,
SURFACES AND CLUSTERS

by

Xinhua Yu

The topics concerning the healing process of metal surfaces; morphological changes of
metal films on insulator substrates; and small metal clusters, have been investigated using
computer simulation techniques.

Surface diffusion is usually the dominant dynamic process by which a solid surface
reaches equilibrium, particularly for metals at temperatures below melting. Above the
roughening temperature, the surface profile and surface tension are continuous functions
of surface geometry. At equilibrium, the chemical potential which is related to the surface
curvature, must be equal everywhere on the surface. Below the roughening temperature,
the surface tension has cusp points that have to be propertly treated in order to derive the
dynamic equations. Two different approaches are utilized and the corresponding dynamic
equations result in the formation of facets on an initially sinusoidal surface profile.

During the formation of thin metal films on an insulator substrate, a sequence of
morphologies occur as deposition continues. An interrupted coalescence model(ICM)
provides a good explanation of the high percolation coverage observed in experiments.
The interruption of coalescence is explained as the combined consequence of

inhomogeneous substrate pinning and thermal fluctuations.

Gold clusters of size N>lO, undergo a fairly sharp solid-liquid phase transition. The
transition temperature decreases as the cluster size decreases. When N is sufficiently
Iarge(>200), the surface atoms have diffirsive behavior in the solid phase. We use the
radius of gyration as a shape parameter to study the shape change of clusters. The radius

of gyration behaves distinctly in various size and temperature regions.

DEDICATION

For those long working days and even longer weekends
while working on the dissertation, with love, I dedicate this
work to Zhixia and Albert.

iv

ACKNOWLEDGMENTS

First of all, I would like to express my gratitude to my advisor, Professor Philip M.
Duxbury whose insightful guidance and broad knowledge in condensed matter physics
help this work move forward smoothly. Without his steady support and encouragement,
this work would never have been initiated. I am also benefited from numerous valuable
lectures and seminars given by Professors M.Thorpe, S.D.Mahanti, David Tomanek, and
T.Kaplan. I sincerely thank Professors M.A.Dubson, B.Sherrill, D.Stump who serve on
my committee.

Here are a few of many families and friends I like to thank: My parents whose love and
support inspire me to do more things; My host family, Wally and Marilyn Baird who
generously accept me as part of their family; Caihua Gao, my mother-in-Iaw who came
from China to help baby-sit Albert; Tie Lan and Ying Liang, Jiahwa Yeh, Yiliang Liu and
Vincent Ma with whom my family did many memorable wonderful activities together;
Qing Yang and Fei Deng whose friendship also enriches my life at MSU.

The financial supports of the Department of Energy and the Center for Fundamental

Material Research are gratefully acknowledged.

Tables of Content page

LIST OF TABLES

LIST OF FIGURES

CHAPTER ONE:

CHAPTER TWO:

CHAPTER THREE:

CHAPTER FOUR:

O....0.....IOOOOOOOOOOOOCOOOOOOOO ...... OOViii

..OOOOOOOOOOOOOOOOOO. OOOOOOOOOOOOOOOOOOOOO 1x

INTRODUCTION 0 O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O 1
1.1 Complexity of Surfaces.. ...... . ......... 1
1 O 2 overViews O O O O O O O O O O O O O O O O O I O O O 0000000000 3

THE ALTERATION OF A SURFACE PROFILE

DUE TO SURFACE DIFFUSION........... ........ 11
2.1 Introduction................... ........ 11
2.2 Surface Diffusion Equation... ........ ..17

2.3 Surface Diffusion Below
the Roughening Temperature... .......... 30
2.3.1 The Anisotropic Effect of
Surface Tension........... ....... 30
2.3.2 "Miscut" Model........... ...... ..33
MORPHOLOGY OF THIN METAL FILMS.............44
3.1 Introduction.................... ....... 44
3.2 Morphology of Thin Metal Films. ....... .45
3.3 The Effect of Inhomogeneity in
Substrates................ ........ .....63

MOLECULAR DYNAMICS STUDY OF

GOLD CLUSTERS........... ................... 77
4.1 Introduction... ........................ 77
4.2 Many-body Interaction in metals. ....... 80
4.3 Solid-liquid Phase Transition... ....... 84

4.4 Kinetics of Non-equilibrium

Gold Clusters................... ...... 106
CHAPTER FIVE: TRANSPORT PROPERTIES ON ZD
PERCOLATION LATTICES......................124
5.1 Introduction............... ......... .124
5.1.1 Anderson localization .......... 124
5.1.2 Scaling theory of

localization.............. ..... 126

5.1.3 Percolation theory.............130

5.2 Calculations of Conductance. ...... ...132
5.2.1 Kubo's Formula......... ..... ...132
5.2.2 Landauer's Formula.... ..... ....135

5.3 Numerical Calculations...............137
5.3.1 Quantum percolation

threshold......................142

5.3.2 Scaling of conductance.........146

5.3.3 Conductance distribution.......152

5.4 Discussions and Open Questions. ...... 156
FINAL REMARKS.................................. ........ ..159
APPENDIX A.. ............................................. 163
REFERENCES....................................... ...... ..168

vfi

LISTOF TABLESOOCOOOOO0.0...OOOOOOOOOOOOOOOOOOOOOOOOOOOOPage
3.1 Critical exponents of ICM calculation and ZD
percolation case................................ ..... 62

4.1 Unit conversion in Molecular Dynamics simulation ..... 93

vfii

LIST OF FIGURES.......................... ..... .. ........ page
2.1 Wulff construction near a cusp point.................13
2.2 A perturbed spherical surface................ ...... ..20
2.3 Two touched droplets............................ ..... 22
2.4 Geometrical illustration of two droplets

with a connecting neck.......................... ..... 23

2.5 Definition of w and 9

......... ......... . ........... ..25

2.6 Three coalescence stages.... ..................... ....27
2.7 Temporal growth of neck size y..... .................. 28
2.8 Surface energy v.s. droplet size.... ................. 29

2.9 Surface profile through solving
equation (2.29)....... ............................... 34

2.10 Temporal growth of facets............................35

2.11 Geometrical picture of miscut model .................. 37
2.12 Numerical solution to equation (2.38)...... ....... ...41
2.13 Temporal changes of facet size and position.. ........ 42

3.1 SEM photographs of indium evaporated onto 8102.......47
3.2 The droplet configuration of MFM

from simulation ................ .................... .52
3.3 Droplet size distribution for MFM. ....... . ...... .....53
3.4 The droplets configurations

generated from the ICM(D=3)........ ........ . ....... ..55
3.5 The percolation coverage Pc v.s. Rc/R0(D=3). ......... 56

3.6 The droplet configuration from the ICM(D=2) .......... 58

The percolation coverage v.s. Rc/Ro (D=2) ............ 59

The cluster size distribution

near percolation point..................... ..... .....61
The experimental results of the Rc dependence

on substrate temperature..................... ........ 66
The geometry of two droplets with

an interior boundary...................... ........... 68
The surface and interface energy

v.s. droplet size....................................71
A typical caloric curve(N=13)............. ........... 86
A schemetic figure showing how the free energy

changes with temperature and floppiness .............. 89
Experimental measurement of gold cluster melting
temperatures as a function of sizes ................ 91
(a)The caloric curves of clusters N=55, 177, 381

and 725; (b)The specific heat v.s. T........... ...... 95
The size dependence of melting temperatures and
comparison with experiment................. ..... .....96
The mean square displacement (6d)2 as a function of
simulation time for cluster N3725...... .............. 98
Diffusion constant D as a function of

T in cluster N=725............. ..................... 100
The mean square displacements of surface

and bulk atoms at various temperatures .............. 102

5.1

The surface diffusion constant D, as a function

of temperature in cluster N=725...... ............... 103
The mean square displacement behavior at various
temperatures in cluster N=381....... ................ 105
A typical plot for initial and final

configurations of cluster N=41 .....................108
rg as a function of time................... ......... 110
Distribution of completion time tc at T=300K

in cluster N=32.....................................111
The behaviors of rg(t) at various temperatures

and system sizes....................................114
The decaying behavior of rg above

the melting temperature........... .................. 119
The decay constant I v.s.temperature T .............. 120
Phase diagram showing distinct

decay behavior of rg ...................... ......... 123
A 32x32 percolation square lattice embedded

in a perfect square lattice belt of width 32

and infinite length..................... ...... . ..... 140
Conductance g as a function of system size.L ........ 144
g as a function of system size in the case

of site percolation.................. ............... 145
A log-linear plot of g as a function of

p at fixed system size..............................147

Linear-log plot of g v.s. L at small disorder ..... ..150

log(g) as a function of L in the strong

localization region........................... ...... 151
The distribution of g in the strong

localization limit............................ ...... 154
The distribution of g in the medium

disorder region.....................................155
Log-linear plot of g as a function of L

at p=005 O. 0000000000000 ...OOOOOOOOOOOOOOOOOO ....... 158

xfi

CHAPTER ONE

INTRODUCTION

W

During the past decade, surface research has clearly
shifted its interest from the macroscopic to the microscopic
scale, mainly thanks to the innovation and widening
applications of electron microscopes; a wealth of novel
experimental techniques and theoretical methods have been
applied and developed successfully. Many subjects such as
catalysis; microelectronic devices and contacts;
lubrication; resistance to erosion; creep and intragranular
fracture etc.; are brought alive again after years of
inefficiency and frustrations. To understand those complex
problems, one must address several fundamental problems:
atomic and electronic structure; diffusion along surfaces or
across interfaces; the energy and chemistry of surface
regions and the response to external forces. Existing
techniques for investigating surfaces have reached maturity
and are increasingly being applied to systems of practical
relevance. On the experimental side, a variety of
microscopes are providing structural and electronic
information about surfaces at a number of complementary

length scales. 0n the theoretical side, phenomenological

2

theories for total energy calculations have had much success
in solving very complex structural problems on the atomic
scale. More sophisticated theories have been introduced to
direct calculations of motions of surface particles by
comparison with experimental data.

To gain comprehension of surface phenomena, it is
critical to understand the complexity aroused from the
termination of a crystal by a surface. The manifestation of
the complexity can be studied from two different aspects:
One is to study the atomic arrangement in the surface
region. Two particularly interesting effects may take place
due to the presence of a surface. (1)the distance between
the surface plane and those parallel to the surface deviate
from their bulk values; (2) some lateral structural change
may also occur to the atoms on certain surfaces(for
instance, Au:110 surface) and this phenomena is called
surface reconstruction which has caught much attention
lately. Another aspect is dynamic properties, for instance,
vibrations and diffusions, which are of considerable
importance to the elucidation of surface phenomena and rich
experimental results. The dynamic problems of surfaces are
usually very complicated and not analytically solvable
because of nonlinearity in dynamic equations and thus a
certain approximation like the harmonic approximation is

usually made. The surface atoms are less bonded than their

bulk counterparts and therefore their mean square amplitudes
are significantly greater than in the bulk. Even at very low
temperature the surfaces appear to be disordered while atoms
in the bulk are very much in order. Often, one needs to
study the migration of particles on a surface; for example
in wetting phenomena, surface flattening and roughening and
thin film formation. Other approaches have to be used to
tackle these difficult problems. Because many of those
problems are associated with phase transitions and critical
phenomena, scaling and renormalization are powerful tools in
the investigation. Computer simulation is adding a new
dimension to scientific investigation, establishing a role
of equal importance with the traditional approaches of
theory and experiment. As the speed and accuracy of
arithmetic computers improve tremendously, it is feasible
and advantageous to carry out quantitative investigation
through detailed computer simulation. It is thus no surprise
that many portions of this work are numerical analyses or
computer simulations.

In the next section, I will give a brief introduction to
the contents of each chapter of the thesis and some key

concepts are outlined.

The first part of the thesis studies the morphology of
thin metal films, surfaces and clusters. The equilibrium
geometry of a finite crystal is determined by many factors
among which the crystal structure, temperature and crystal
size are considered the most important. At any temperature,
the shape of a finite crystal can be found theoretically by
the Wulff construction( Wulff, 1901) provided that the
surface tension 7 as a function of crystalline orientation 9
is known. At finite temperature, the total free energy of
the system reaches a minimum at its equilibrium shape. If
the temperature is sufficiently low, it is expected that the
geometry does not change much from its ground state except
at those positions like corners and edges where atoms are
easily activated. As temperature increases, the surface
atoms start to diffuse prior to the bulk ones because of
their reduced binding. One fundamental question is how those
diffusing surface atoms alter the non-equilibrium shape of
crystals. This question is to be discussed in chapter two.
One approach is to view the surface to be a continuous
function f(x,y,z)=0. The diffusion of surface particles is
caused by the gradient of chemical potential u which is
proportional to the curvature K.of the surface if the
surface tension 7 is isotropic. As a result, a differential
equation describing the process can be derived. The surface

diffusion equation has many successful applications in

5

studying systems at temperature T~0.75Tm, the bulk melting
temperature. Nevertheless, as a matter of fact, 7 is
strongly dependent on surface orientation at low
temperatures and particularly has a cusp singularity along
crystalline planes of high symmetry:y(€)=ro+7,|fl. To
overcome this difficulty, two methods are introduced and
alternative equations are derived. One method originally
developed by Bonzel etc.(Bonzel, 1984,1986) is to de-
singularize the cusp point by replacing the small angle 9
with ‘/92+6(2, -60 where 90 is a small constant. As a result,
the second derivative of 7 exists. The other method assumes
that the surface profile has a miscut angle with the high
symmetry lattice plane. The alteration of surface shape is
due to the motion of steps. Lancon and Villain(Lancon and
Villain, 1990) studied the chemical potential of steps and
obtained a different version of surface dynamics. Because of
the nonlinearity of the equations governing those dynamic
processes, their solutions are attained via finite
difference methods.

The work on the formation of thin metal films( Chapter
Three) on an insulating substrate was motivated by the
experiments by Dubson and Jeffers( Dubson etc. 1994). When
metal vapor condenses on a clean non-wetting insulating
substrate to form a thin film, a sequence of morphological

changes occurs as the thickness of the film increases, going

6

from isolated compact islands, then a "wormy" structure, to
a percolating structure and finally a continuous thin film.
It is observed that an abnormally high percolating coverage
often occurs. Based upon the complete coalescence model by
Meakin and Family(MFM) ( Family and Meakin, 1988,1989), we
proposed the interrupted coalescence model(Yu etc.,1991), in
which two droplets on the substrate do not coalesce as long
as the sizes of both touched droplets reach a critical
dimension RC. Rc depends on deposition rate, substrate
temperature, as well as the substrate structure. Because of
the presence of RC, this "interrupted coalescence" model
goes through stages of early complete coalescence, crossover
to wormy structure, percolation and hole filling as in
experiments. Our computer simulation mimicking the
experiment showed that the coverage at the percolation point
increases with increasing RC. We also calculated the cluster
size distribution and the critical indices near percolation
and found that this model is in the same universality class
as the conventional percolation model(Satuffer,1985).
Furthermore, it is shown that the inhomogeneity of the
substrate may be the mechanism leading to the interruption
of coalescence. This is illustrated by a calculation of the
total energy of two droplets with a grain boundary in their

interior.

In Chapter Four, I deal with smaller physical systems,
metal clusters, using the molecular dynamics simulation
method that is becoming an increasingly powerful tool in
scientific investigation. The simple idea of molecular
dynamics is that the detailed information such as position,
velocity etc. of each atom in a system is calculated as a
function of time assuming that each atom complies with
Newtonian dynamics. From this information, thermal dynamics
quantities can be calculated straightforwardly. In a metal,
many experiments reveal that atoms have many-body
interactions, therefore it is critical to have a good
potential to describe the metal. We have chosen an embedded-
atom potential(Gupta, 1981) which gives good predictions for
many bulk properties of transition and noble metals. Because
of the outstanding surface properties of gold and the
availability of experimental data, we selected gold clusters
for study. The calculations of temperature and total energy
show that the liquid-solid phase transition is clearly
identified in clusters with size greater than about 10
atoms. The melting temperature Tm decreases as system size
decreases and obeys the following semi-empirical

relationship(Buffat,1976;Borel,1981):

TAR): Tm(oo)(l—%) (1.1)

8

where Rc is constant and determined by the bulk properties
of gold. By calculating the mean square deviation of both
surface and bulk atoms, we found that the surface atoms have
diffusive behavior at temperature below melting for large
clusters(N>200) in comparison with the absence of surface
diffusion for small size clusters( N<150). The dynamics is
studied by looking into the shape of non-equilibrium gold
clusters and the radius of gyration rg represents the shape
parameter. We observe how r9 decays from its initial value.
Above the melting temperature, the correlation between atoms
is weak and r9 decays roughly exponentially with large
fluctuations. Below Tm, the decay of rg depends on the

cluster size N. When N is small (less than 150), r reaches

9
a meta-stable value quickly with typical time te and remains
almost unchanged for a long time tc and suddenly decreases
to a smaller value within a very short time tq. The three
time scales have the following relations: teth<<tc. When N
is large ( greater than 250), and the temperature is close
to Tm, rg decays like smaller clusters. But when T is
considerably away from the Tm, rg decays algebraically
because of the dominance of surface diffusion in large
clusters.

Chapter Five of the thesis describes work on quantum

transport of two-dimensional percolating lattices. Electrons

are permitted to hop between two bonded sites on a

9

percolating two dimensional lattice. Therefore the
Hamiltonian describing the system is simply written

as:(Anderson,P.W., 1958)
H=Zt0(ara: +ar+aj)+zsiar+ai (1'2)
1'.) i
where i,j are nearest neighbor lattice sites and szl for a

present bond and zero for an absent one and 61 is the site
energy. From percolation theory, when the bond probability p
is less than a critical value pc( for a 2d square lattice
pc=0.5) there exists no infinitely large clusters. For a
system of large but finite dimension, the largest cluster
size has a size of the same order as the system. If all
present bonds are good conductors and a voltage is applied
to two sides of the system, a current will flow when p>pc,
and no current flows for p<pc. Because of quantum effects,
an electron is more likely to scatter backward and more
transmitting bonds are required for an electron to transmit
across the sample. In other word, the quantum threshold p9
below which no current flows through the system ought to be
larger than pc. To find p9 and investigate the scaling of
conductance, we employed two methods in the calculation of
conductance: i) the linear response theory and ii) transfer
matrix and Landauer formula. In the linear response theory,
we introduce the Green's function G¥=(E-Hiin)'1 where E is

the Fermi energy and n is a small positive number. The

10

conductance can be expressed as the product of two Green's
functions which are expressed in matrix form and are easily
calculated. The conductance Gc is related to electron
scattering by disorder via Landauer's formula:
Gc=(e2/h)tr(tt*) where t is the transfer matrix of the
disordered region. Our calculations show that p9 is one for
2d percolation lattices. This result supports the prediction
of weak localization that any finite disorder will cause
localization of the electron wavefunction in two dimensional
systems.

Finally I summarize the work in the thesis and discuss
several interesting future projects stimulated by the

analyses described here.

CHAPTER TWO

THE ALTERATION OF A SURFACE PROFILE
DUE TO SURFACE DIFFUSION

gamma

An isolated system, say a heated metal crystal, tends to
alter its shape in order to minimize its total surface free
energy. The common processes by which the shape change is
realized are evaporation-condensation, volume diffusion, and
surface diffusion. The process of evaporation-condensation
(Mullins, 1957, 1959) occurs when the chemical potential of
the system u differs from that of its co-existing vapor Ho.
Evaporation- condensation does not conserve the total number
of atoms in the system which gains atoms when its chemical
potential is lower and loses atoms when its chemical
potential is higher than the vapor's. At finite temperature,
there exist bulk vacancies and defects which allow the bulk
atoms to move, in which case volume diffusion is important.
The ratio of vacancy number to total atoms is approximated
to be exp(-Eb/kT) where Eb denotes the approximate energy
needed to create a vacancy. For most metals, £5 is of order
several electron volts, i.e. 104-105K and thus the ratio is
extremely small at room temperature. Surface diffusion
involves surface particles which move along the surface of

the material. Due to the reduced bonding at the surfaces,

12

diffusion processes are much more probable at surfaces than
in the bulk. For a small metal particle with size 0.1 to 10
um at temperature of about 0.75Th (T5 is the bulk melting
temperature) the surface diffusion is considered to be a
dominant factor governing the morphological evolution.

The equilibrium properties of a surface can be described
by a few thermodynamic quantities like surface tension 7
defined as the surface free energy per unit area. Generally
y is anisotropic and depends on the orientation of the
surface. The lowest energy profile of a surface is then non-

trivial and the total surface free energy of a finite system

fizfiy-dA (2.1)
A

has to be minimized under the constraint that the volume is
constant. The equilibrium shape of a crystal of fixed volume
can be determined from the Wulff construction (Wulff,1901;
Herring,1951) which utilizes the result that the surface
tension 7 of any face is proportional to the distance of
that face from the center of the equilibrium shape.
Therefore, the process is described as follows: give a polar
plot of surface tension 7(6); draw a ray from the center at
any direction 6 and the line intersects with the y(9)

surface at a point; construct a plane perpendicular to the
ray(Wulff plane); the crystal shape is the envelope of all

those planes. One particularly

Q!

-P'

13

 

 

Figure2.l:1hepreeeneeofacusppointin1(9)willleadtothe
facets on the equilibrium shape through Wulfi‘eonstrucuon whose
procedure was briefly described in the text. The polar plot near 9
=0 is forthe function 7(9)-yo+1,|0|.

14

interesting consequence of the Wulff construction is that if
7(9) has cusp point, a finite flat area (facet) will be
produced in that direction. £1gn;g_zL1 shows why this is the
case. From this figure, it can been seen that the envelope
formed by plane SP is the plane P'P as we move S along O'Q.
The surface atoms are highly correlated at very low
temperature. As the temperature increases, a two dimensional
surface undergoes a roughening phase transition(Week,1980).
Let h(x,y) denote the height of a surface with respect to a
proper reference plane (e.g. the perfect surface at T=0),
where (x,y) labels the position on the plane. The roughness
of a surface is defined as the spatial correlation function:
(w = (Rm-(hf) (2.2)
Where <> stands for thermal average. Without loosing
generality, it is assumed that above the roughening
temperature TR, the leading terms in a Hamiltonian
describing the surface are(Fisher,M.E., 1988)
H=1/2j¢cdy[J(Vh)2+gh2] (2.3)
when both Vh and h are small. Here J and g are constants
relating to surface tension and external force respectively.
Then (Ah? is calculated according to a Boltzmann
distribution, giving (AhfcrhrL, where L is the surface
dimension. Below TR, the Hamiltonian ought to include terms
which favor h being integers. Using either solid-on-

solid(S.O.S) or the sine-Gordon model, it is theoretically

15

proved that the roughness of a surface is only a function of
temperature and independent of system size at low
temperatures. The surface free energy at T-TR has a cusp
singularity point along directions of high symmetry.

The experimental methods(Schommers etc.,1987) used in
studying surface diffusion are roughly divided into three
groups: a) direct observation of diffusing atoms using field
ion microscopy (FIM) which records the motion of individual
atoms or clusters on a small metal surface area and electron
microscopy which observes the migration of surface atoms and
simultaneously the detailed surface structure. The diffusion
constant D can be directly obtained from FIM data using
<rflj)=4lh) and the temperature dependence of D gives the
magnitude of barriers over which the diffusion process is
thermally activated, i.e.

D=Doexp(-E,/kT) (2.4)
where E8 is a surface activation energy. b) Concentration
gradient methods such as scanning AES(Auger electron
spectroscopy) analysis that studies hetero-diffusion and
laser light diffraction which attains information on self-
diffusion on metal surfaces. The laser light diffraction
experiments measure the amplitude of a surface profile at
constant temperature as a function of annealing time by
recording the intensity distribution of laser-light

diffraction generated by the periodic profile. c)

16

Equilibrium methods which include the field emission
fluctuation method, neutron scattering and He-atom
scattering, chiefly studying the diffusion of adsorbates.
These experiments employ inelastic scattering to measure
space and the time dependent self-correlation function
C(r,t) or its Fourier transform S(q,m).

Theoretical study of surface diffusion relies on the
continuum model which assumes a system surface is described
as a continuous function z(r,t). The model proposed that the
surface particles diffuse because the chemical potential u
varies along the surface. The mechanism is very similar to
the concept of electric charges moving along a conducting
surface when a conductor is not in equilibrium (equi-
potential condition). It will be seen that the continuum
model describes self-diffusion above TR quite successfully,
but has serious problems when the temperature is below TR.
More detailed study is possible with either molecular
dynamics(MD) or Monte Carlo(MC) simulation techniques, which
requires a suitable atom-atom interaction.

In this chapter, I will introduce previous work on
surface diffusion relating to our later discussions and
present some numerical calculation results. In the following
section, a diffusion-like equation describing the dynamics
of surface profile will be derived using a continuum model

and some simple applications will be presented. In section

17

2‘; I will discuss surface diffusion below the roughening
temperature TR. We discuss the anisotropic model first and
then the miscut model. The conclusions are summarized also

in this section.

WW
Let us start with the expression for the free energy of

a bulk-surface system with surface profile z(x,y). The
expression is written in the form:
F=pollm(x.y)+fiabcdya(z;.z;) (2.5)
Here the first term is the bulk free energy and the second
one is the surface free energy and 0 relates to surface

tension through

a: y(z',,z'y) 1+:r'2 +z'2 (2.6)

3 y

The chemical potential is obtained from p=QéF/62,
yielding:

Wong flips) (2.7)

a: 3: ay 52;

 

If y is constant, one has
a z' 6 23

fl'.“ ="Q}’ ———L—"—+—' ', (2-3)
0 6% /]+z':+z': @Jl+z';+z'i

It has been shown that the term in the parentheses on the

 

right hand side is exactly equal to -(K1+K2)=-K, here K1
and K2 are the two principal curvatures of the surface.

Therefore we have the Gibbs-Thompson relationship:

18

fl-uo=QIK (2.9)

For better illustration, we consider one-dimensional
isotropic surfaces of height z(x,t) and derive the dynamic
equation first, then extend to anisotropic and higher
dimensional systems. Suppose z(x,t) coexists with its vapor
of chemical potential “0 and that the bulk chemical
potential is also “0 to suppress the evaporation-
condensation process and surface diffusion is the dominant
dynamic process to alter the surface profile. From equation
(2.9), if the curvature varies on the surface, the surface
particles undergo a "field", the gradient of chemical
potential which drives them to drift. The drift velocity
ought to be proportional to the magnitude of the field,
reading

v =-M —%2K—

V = 2.10
. . .# k,T a ( )

Here the Einstein relation Adszlg/kBT has been used and M8
and Da denote the mobility and surface diffusion constant
respectively. The corresponding surface current is the

product of drift velocity and area density v, reading:

L>fl2véK.
‘/=_ §_____ 2.11
kBT & ( )

Applying the continuity equation to the surface,
£31,913

:0 .
as a (2 12)

19

where p denotes the particle number variation per unit
surface area, one easily obtains the distance change along

its normal by multiplying atomic volume (2
a _ a _ 1),;ozv52K
— a 1:32“ a2

 

(2.13)

It is easily generalized to higher dimensions by replacing

0'2ch%2 with VfK, which gives

%=BV§K (2.14)

where B=D37flzvlkBT and K is sum of two principal curvatures
at a surface point. When 7 is not isotropic and dependent on
surface orientation 9, the modified Gibbs-Thompson relation

(Herring,1952)is:
(77 ¢?y
p-po=fl(7+—2-]K,+Q[7+-—2-)K2 (2.15)

I 2
on the condition that the second derivative of 7 exists.
Here 61 and 62 are angles between the normal at a surface
point and two principal planes. Correspondingly the surface

diffusion equation becomes:
%=§V‘2[[7+%)K,+[y+%)l(z] (2-15)
The equilibrium shape when 7 is isotropic is trivial as
seen by setting n to zero, and it is simple to show that
the equilibrium shape is a circle in two dimensions and a
sphere in three dimensions for a closed surface. Now we add

a perturbative term to a two dimensional equilibrium surface

and calculate how the shape relaxes under the perturbation.

 

20

 

l’l .
/ 6R \
\
’ l
l
l
l R 1
\
/
\ /
\
~ ”

Figure 2.2: A perturbed spherical surface. The perturbation
magnitude is much smaller than the radius of the sphere so that
one can treat the problem as small perturbation.

21

Adding a small circumferential perturbation R0A(6,¢) to the
equilibrium shape r=Ro with A(0,¢) <<l(see W) . The
surface Laplacian is written in terms of spherical

coordinates (r,6,¢) to be(Brailsford and Gjostein, 1975):

 

1 a a 1 a?
V2=—sin —sin9—+ — 2.17
’ R0 39 59 sinzaapz] ( )
A(9,¢) can be expanded:
#242,449») (2.18)
1,»:

where Yim is the standard spherical harmonic function. To

the first order of A, the curvature is also written as:

K=KO+7I:;Z(I+2)(I-1)A,MY,M (2.19)

Ln:
It is immediate to obtain the differential equation for each

A m component

Alm=-%l(l+2XIz'-1)Alm (2.20)

whose solution is

Aim“) = A[m(O)CXP(—I/ Tim) (2 ° 21)

where 7':=-I%I(I+2)(12-1). We have assumed that at equilibrium
0

A,m(oo)=0. The result shows that the relaxation time 11m for
mode {1,m) in proportional with R3. When l>0, thn-(Ro/l)4,
which is approximately the fourth power of the wavelength.

Therefore, short wave perturbations decay much faster than

long wave ones.

22

 

 

(a)

(b)

 

 

 

 

 

Figure 2.3: Two touched droplets as an initial profile for
numerical simulations. The neck area has been smoothened using

a small circle of radius r. (a): overall profile (b): magnification of
neck area.

23

 

 

Figure 2.4: Geometrical illustration of two droplets with a
eonnectingneckshowingtherelationshipbetweenrandy:
r-y3/2(l-y)~y’/2

24

One of the very interesting applications of equation
(2.14) is the estimation of the complete coalescence time of
two droplets of size R with a very small neck connecting
them (£1gn;g_z‘1). The neck area is approximately the shape
of a circle of radius r. The geometrical relationship
between y, the neck size and r is better illustrated in
£1gg;g_2‘4. Actually, the growth in the neck area is
extremely rapid when y is small because of the large
curvature. Here we roughly estimate how fast y grows when r
is small. Using R=1, we can write down an approximate
dynamic equation for y: y=B§K/&2. Notice that r5y2/2 and
3K/o‘szz(2/r)/(n2r2), for very small y, therefore,
y=B(l6/Ir2)y" or tzO.ly7/B. Suppose y=0.2, tz1.2x10’5/B. In
order to perform numerical calculations, one needs to
convert the equation (2.13) into dimensionless form, which
is given by:

—=fo (2.22)

where N=n/R0,R=r/R0,S=s/RO,K=KRO,r=Bt/R:. It is convenient
to express K and the Laplacian in terms of polar

coordinates (R,8).

(23 41:2,)?2 (99.359 '

From the illustration diagram Eiggre 2.5, the curvature is
defined as dV/ék. The relationship between 8 and w is

w=0+7r/2—tan"(R'/R) (2.24)

25

 

 

 

Figure 2.5: Illustration of a surface segrnent. showing the
definitiooofvandeandtheirrelationship,fmmwhiehdie
curvatureeanbeexpreuedinpolareoordinateflrfi).

 

26

Therefore, the equation upon which numerical calculation is

performed can be written as:

at 15 1 a:
__.___ _ 2.25
(9? R50 [8+sz ( ’

where R=dQ/69. In same notation, x is expressed as
x_ R2 +21?2 —RR”
(R2 +R‘2 )3/2

The surface shape is described by the function R(9). I

 

(2.26)

employed the standard finite difference method to evaluate
the derivatives. Convergence requires that the von Neuman
stability criteria is

61am“ (2.27)
£1gg;g_2&§ shows how the surface shape evolves with time.
Because of the extremely large curvature near the neck area
in the early stage, the neck increments very rapidly. The
plot shows three stages of shape change: early neck growth,
neck elimination( at time t=0.17) and complete coalescence(
at time scale 1.0). As a matter of fact, the neck increments
from an initial y=0.05 to 0.6 within time scale 10’2. The
process slows down significantly when y20.7. Figure 2.7
explicitly demonstrates the temporal growth of neck, showing
that the growth speed is very large at the beginning and
slows down, eventually turning into an equilibrium shape.
Our calculation results are consistent with experimental
observations(Nichols and Mullins,1965). During the process

of coalscence of two droplets, the surface energy decreases

CO
Q) D

 

Figure 2.6. The droplets shape as a function of simulation time.
The figure shows three different coalescence stages: (a): early
stage: 3:055, t=0.008; (b): neck elimination: )=O.90, t=0.17 and
(c): full coalscence? }'=1.42, t=0.90.

alze y

1.5

1.0

28

 

 

 

 

 

 

I I I 1 I
)- d
l— —
- -(
- a
-— '1
l I l l J l J l l l l I l I l l J l l l I l l l L L
0.00 0.25 0.00 0.75 I.” 1.20

timet

Figure2.7:Theternporalgrowthoftheneeksizey.'l'hegrowthis
Wynpidintheaflysmgeanditslowsdownafier
rachingaeertainsinmarticularlyafierneckelinunatiou.

11.0

10.6

10.0

9.0

10.5

10.0

29

 

TIrTlUVIU'TTrrlIIrrIITUU

.
d
‘
.—
q
q
.
d
—
d
d
d
.
—

(a): r,-1.0.r.-1.0

IlLlllLLlllLllllllllllll

 

 

 

 

 

Figure 2.8: Surface energy as a function of droplet sizes. The
surface energy is directly related to total area or length of the
droplet surfaces. The droplet size is represented by the distance
between two farthest point on the surface. (a): two equal size
droplets (b): two droplets having different sizes.

: : : : l : s : . ' : : . : J : J.
3: ile ila zlo :
:— (b): r,-o.75.r.-1.2 -:
' i
C ..
C I
f' "I.
b l l l l l l l j l I 1 l l 4 l L 1‘
.4 16 1.5 20

30

as the system approaches equilibrium. We studied how the
surface energy is altered as the system size and the
distance distance between two droplets various. This is
shown in Eigg;g_zL§. At the early stage, the relaxation
process in the neck area is so rapid that the system size
does not change at all while surface energy decreases

rapidly.

WW
2.3.1. Anisotropic effect of surface tension

The continuum surface diffusion theory has had much
success especially in the temperature region T>0.75Th close
to the melting temperature. As T is decreased to below TR,
the crystal tends to develop facets on its surface which has
been intensively studied by experimental physicists during
the last decade(Jayaprakash etc.,1983). Most experiments
were done on metals at room temperature, much lower than
their melting temperature. The samples are prepared with
periodic surface profiles and exposed directly to an
electron microscope which records the shape as a function of
time. It is observed that facets start to take shape at tops
and bottoms of the initial curve. At the beginning the
process is quite fast but slows down as the size of facets
becomes large. From the isotropic continuum model, this

phenomena is not expected to happen because a sine curve is

31

a eigenfunction of the equation so that the sinusoidal shape
is preserved and only its amplitude decays exponentially.

As pointed out by Landau (Landau, 1986), the surface
tension 7 of a planar solid surface depends on the
crystallographic orientation 0 of the sample. Thus as we
calculate chemical potential, the modified Gibbs-Thompson
relationship is utilized. When 7V<7k, the surface tension 7
has a cusp point at 9:0. Near this point:

7(€)=7o+r.l6i (2.28)

When substituting equation (2.28) into equation (2.15),
immediately we have the difficulty that the second
derivative of 1(9) is singular(a 5-function) at 9=0. The
dilemma was first circumvented by Bonzel and coworkers
(Bonzel etc.1984; Preuss,1986) mathematically by
regularizing the 6-function.

For better illustration, again we consider a one—
dimensional surface z(x). The equation of motion for the

surface is written as:

-——[7(6’)+

 

a?

86;]Iqx) (2.29)
The above equation is greatly simplified using small slope
approximation z'<<1 and reads

. B
z=-—17"(x) (2.30)

0

where prime represents differentiation with respect to x and

32
'Kx)=[r(9)+-§0,Z]Z" (2.31)

Here 9 is also a function of x and Osz' for small slope. To
regularize the singularity of the second derivative of y,
the key substitution is to replace 9 with JON-6% —00, with
parameter 90 being very small. This substitution does not
change the:y(9) much when 9>>90 but provides finite second

derivative of 7(6) near 680. Equation (2. 31) turns into:
79;
17(x)=[70+7,(‘/z'2 +93- 60 )+(z —-1—0¢:—-—£)3,2]z" (2.32)

by replacing y with 70+7,(W-(;o). To analyze the overall
behavior of this equation, we regard z', z", z'" and 2"" as
small so that the second term on the right hand side is
always small compared with the first term. Combining
equation (2.30) with equation (2.32), we have the simplified

differential equation:

 

(2.33)

The first term is the classical surface healing term
obtained by Mullins(Mullins,1957) and the second term
embodies the effect of anisotropic surface tension and
vanishes if Qy=0 and z'to. This equation can be discussed
in two specific regions. First in the region where z'SGO,
equation (2.33) becomes

. B
zz-BZION___)/Lzllll (2.34)

6070

33

Typically 75/n320J, while 90 is taken to be less than 10’3,
therefore, the second term on the right hand side dominates.
In the region that z'>>90, the second term is negligible
compared to the first term. The speed of decaying in the
regions of z'SBo(flat region) is much faster than in the
other region, therefore an appreciable shape change of the
surface profile occurs at the tops and bottoms where facets
are initially developed.

The direct solution of equation(2.29) through numerical
methods is obtained by imposing periodic boundary conditions
and a sinusoidal form as an initial condition. Eigg;g_z;2
presents a sequence of surface shapes versus time. It is
clear that facets show up at the bottoms and tops of the
curve almost immediately after the simulation is started.
The size of those facets increases with time while the
height decreases and Eigg;g_z‘19 exhibits that time
evolution. Our calculation data shows that the facet width
increases with time t very rapidly at early stages and slows
down after reaching size of order one. Just prior to the
disappearance of facets(due to completion of flattening),
the widening speed increases again. The result is very
similar to the calculation by Lancon and Villain (Lancon and
Villain,1990) who utilized the evaporation-condensation

dynamics and miscut model that is the topic of next

subsection.

34

 

 

 

 

 

 

 

 

 

 

. I I I U l U U I I I I I I I l q
0.10 — _
0.05 — .:
L d
0.00 ‘
' J
-0.0s — L
- 1
-0.10 — _:
l l l l i I l l l l I l I I 1 I

o 2 4 o

I

Figure 2.9: Numerical solution to equation (2.29), showing the
surface profile decaying. The initial shape is sinusoidal. The
formation offaoets occurs at the tops and bottoms while the shape
between facets approximately maintains a sinusoidal solution.

 

width 1r

height h

 

35

 

III

1.6

1.0

[IIII'IIIIIIIIII

0.10

0.00

0.00

0.04

0.02

IIIIIIIIIIIIII'IIIIIIIIII

 

I III I I I I III II IIII'I III I I III III

I
rrsfrlr
N

L

(a)

   

lLllllllLlllllllllllll

 

 

 

000°11L.llllllllllllllllljjlliaa

 

-0.2

0.0 0.2 0.4 0.0 0.0 1.0
time t(arb. unit)

Figure 2.10: The size and height of those facets indicated in
Figure 2.8 as functions of time. (a): the facet size grows with time
rapidly at early stage, slows down after the size reaches the order
of l and vanishes as approaching equilibrium. (b): the vertical
position of facets decreases monotonically.

36

2.3.2 Miscut model:

We start with equation (2.7) to calculate the chemical
potential. When T<TR, the surface consists of well developed
steps and terraces and the surface shape change is due to
the motion of steps. In general, a can be viewed as the sum
of two contributions: 00, the free energy of terraces and
0,, the energy of steps. Further a, can be written as
a, = a, + a, , where 01 is the free energy of non-interacting
steps, which is proportional to the step density; 03 is the
interaction between steps. The existence of 03 is due to the
fact that two adjacent steps tend to be away from each other
such that they could have more space to wander around. A
detailed investigation(Gruber and Mullins, 1967) shows that
03 is proportional to the third power of step density.
Summarizing the above argument, 0 has the following
approximate expression:

a: 00+71|u|+13~y3|u|3 (2.35)
where hd=(ii+f:Y”, representing the step density. If we
substitute equation (2.35) into equation (2.7), the chemical
potential is still singular at u=0. However Lancon and
Villain found that if the lattice planes of high symmetry
are not perfectly aligned with the surface, but instead have
a small "miscut" angle, evaporation-condensation dynamics

leads to the formation of facets on the tops and bottoms of

37

 

 

 

 

Figure 2.11: A geometrical picture of the miscut model, showing
that the lattice planes(x-y plane) and surface(X-Y plane) has a
small angle a. It also shows the steps and terraces.

 

38

an initial sinusoidal surface. Here we apply the model to
the case of surface diffusion dynamics.

In the miscut model, the lattice planes of high symmetry
(e.g. Au 111) are not perfectly parallel to the surface, and
they makes a small angle a with the y-axis of a surface
whose geometric setting is shown in £1gnzg_zL11.
Consequently, each lattice plane will have an intersection
yn at the y-axis. Let yn and z be functions of x only, then
2 and yn have the relationship:

z(x,t)=y,,(x,t)sina+ncosa (2.36)
where n stands for the numbering of lattice planes. The

chemical potential u is calculated by inserting equation

(2.35) into (2.7) and setting 23==a.

ék

and a3 = 3 “((22 +2'2 )3/2 . The equation of

p-p. = -3[§Ia.(z'.a)+a.(zuan] (2.37)

where a, = (2'2 +012)"2

motion for z(x,t) is then written according to surface

diffusion dynamics:

0“: 5
-—=-B——, 2.3
53 ék‘w ( 3)

and here

(024-192)” 7’) 2.2+a2)1/2 "

2 :2
[ a +73 22 +02 ]~" (2.38a)

Again we can analyze equation (2.38) in two distinct
regions. First consider the region where z'sa<<1, the second
term of equation (2.38a) is small compared with the first

term, and hence equation(2.38) is reduced to

39

i=(—B/a)z'"' (2.39)
This is a linear differential equation whose solution is
simply z(t)=z(x,0)exp(-Bk‘t/a) for the initial curve
z(x,0)~sin(kx), and thus the decay is very rapid. When z'>>a,
equation(2.38) turns into:

z°=-B,§,-(|z°(z") (2.40)

The solution to equation (2.40) has been studied by Ozdemir
and Zangwill (Ozdemir and Zangwill, 1990) and it was found
that:

z(x,t)=z(x,0)/(l+}lt) (2.41)
where l is constant. The shape preserving solution can be
fitted very closely to a sinusoidal function. Because of the
extremely large decay speed in the region where 2': 60, the
overall surface profile is broken into two kinds of regions:
one is totally relaxed and located at tOps and bottoms of
the surface(facets); the other regions are those between
facets whose decaying behavior is governed by equation
(2.40).

We studied the numerical solution of equation (2.38)
using finite difference methods. The initial condition has
been set to be sinusoidal again. Because a sine-curve is not
the eigenfunction of equation (2.40), it is anticipated that
the surface profile between two facets will slightly deviate

from its initial shape as seen from Figure 2.12. This plot

40

shows that right after the simulation started, small facets
have been created at the tops and bottoms of the initial
sine curve. The size of facets grows with time while the
vertical position(height) of facets decreases(£ignzg_2&11).
It is noticeable that the behavior described by equation
(2.38) is very similar to that discussed in the previous
subsection. Actually, if we regard the miscut angle a as
the parameter 90, the surface tension in the anisotropic

l/2

model has the form:7=yo+yl(a2+z'2) which will certainly

give the same result as using the substitution
9—) 02+65 —60. The miscut angle (1 thus provides a physical
explanation of 60 introduced in the anisotropic model.

To conclude, a non-equilibrium crystal surface alters
its profile mainly through surface diffusion when the bulk
chemical potential is close to that of the environment the
crystal resides in. Above the roughening temperature, the
diffusion process is governed by a continuum model in which
the chemical potential on the surface is proportional to the
curvature and consequently a diffusion-like equation that
relates the alteration rate and its surface geometry can be
derived to closely describe the diffusion process. When
below the roughening temperature, the diffusion process is
complicated owing to the fact that the surface tension has
cusp singularity points in the direction of the crystal

planes of high symmetry. To circumvent this difficulty, two

0.10

0.05

41

 

I I I I j I T I I T I I j l I I I I I j

 

 

 

 

 

 

 

 

- a
I— —I
L m u
-— —i
- 1
" ‘l
)- III
)- I
i- I
. U d
'- 'l
l l l l l I J 1 J J 1 l l I l l l l l l
.0 0.2 0.4 0 0 0 8
I

Figure 2.12: Numerical solution to equation (2.38), showing the
formation of fleets. We have set a-0.01, 7, / y, = 1.0 and the
initial sinusoidal amplitude is 0.1, wavelengthtlfl

:‘
o

width 1r

height h

42

 

0.5 I I I I I I 1 I I lj r I I l I I I I I I I I I I rj I I

1111

(a)

0.4

0.3

0.2

0.1

IIII'IIIIrIIIIlIIII‘ITII—

 

0(l_0 '
0. 10
0.08
0.00
0.04

0.02

 

 

° l l l l
LlllJllllllLllLllllllllllllu llll llll IILL llll

0.00lllllllLllllLlllllllll 11.1,
-0.02 0.00 0.02 0.04 0.00 0.00 0.10'

time t(arb. unit)

Figure 2.13: (2): Temporal growth of facet width, showing similar
behavior to Figure 2.10. (b): facet height as a function of time, the
decay is monotonic.

43

different models have been presented. The first is to de-
singularize the cusp point by making the second derivative
of the surface tension finite through introducing a small
constant. Consequently, the corresponding differential
equation governing the diffusion process yields solutions
that can explain the experiments that facets are observed at
certain locations. The second model studied the chemical
potential in more detail and proposed that the formation of
facets is possible when the surface is not perfectly
parallel to the lattice plane, but with a small miscut angle
a. If one regards the small parameter 90 as a, the two

models are mathematically very similar.

CHAPTER THREE
MORPHOLOGY OF THIN METAL FILMS

W

From an experimental point of view, thin metal films can
be prepared by a wide variety of techniques. The most common
method of preparing thin films is the so-called evaporation
method(Pashley,196$). This involves condensing the metal
onto a substrate with the vapor source being provided by
evaporation. The first important factor influencing film
growth is the quality of the vacuum in which the evaporation
is carried out. It is quite clear that the various types of
chemically reactive or condensable contaminants which are
present in such a system could have a major influence on the
growth and structure of the films. Modern vacuum technology
has enabled the production of high quality vacuums with
pressures of 10"11 torr(Umbach etc.,1991), which provides a
very clean condition for deposition studies. Many different
techniques are utilized for providing vapor sources, the
simplest being a tungsten or molybdenum spiral which is
heated by passing an electric current through it. The rate
of deposition, and the final film thickness , often needs
careful control, since the structure of a deposited film can

depend considerably upon both parameters. It is common to

45

heat the substrate during deposition. The heating can have
important effects in relation to the cleaning of the
substrate, as well as enhancing recrystallization effects
and general structural changes during growth, which can be
important for increasing grain size. The invention and
improvement of electron microscopy have changed the method
of surface investigation significantly. The electron
microscope has a number of advantages for studies on thin
film growth: (1) the very high magnifications possible,
associated with a resolution below 1nm, allow very detailed
study of localized regions of specimens; (2) selected area
electron diffraction allows detailed correlation between the
image and the electron diffraction pattern; (3) internal
crystallographic structure can be detected and analyzed in
detail; (4) direct observation of the growth of thin films
can be made by carrying out the deposition inside the
microscope. Eiguzg 3.1 gives snapshots of the morphology of
an indium film on a silicon oxide substrate. In the next
section, I will discuss some interesting features of those
diagrams and present a theoretical model to simulate the

formation of these features.

3. Coa esce so a d e co on t etal lms:
As shown from many experiments(Paudit etc.,1982), most

metals do not wet insulating surfaces such as glass,

46

graphite or silicon. In such case a thin film grown by
thermal evaporation passes through a sequence of
morphological changes as the film thickens. When
sufficiently thin, the film consists of isolated, compact
islands which, as more material is deposited, grow and
coalesce into larger, but still compact islands. At some
critical island size, islands that touch no longer fully
coalesce into near-equilibrium compact shapes but instead
form elongated wormlike structures. As growth proceeds,
these ”worms" grow longer and connect to form a percolating
structure, and finally, the channels between worms fill in
to form a continuous, hole free, film. This morphological
sequence is illustrated in Eiggzg&_1‘1 which shows an indium
film evaporated onto a room temperature 5102 substrate. For
this film, the critical area coverage for percolation is
quite high (pc=0.82), a result of the fact that the channels
between worms are very long and narrow compared to the width
of the worms. High pc is a puzzling feature of the
morphology of many metal films grown on warm substrates.

As can be seen in Eiggzg_;‘;, in the early stage of
growth, the film has a distinctive morphology, consisting of
compact islands, all of about the same size, separated by
gaps which are comparable to or smaller than the island
radius. Actually, the gaps are filled with a population of

smaller islands which do not show well in these SEM

47

 

Figure 3.]: SEM images of indium on an $02 substrate. showing
complete coalescence; partial

difierent stages: the early
coalescence; wonny structure; percolation and hole filling.

48

photographs. During film growth, islands are always growing
because of deposition of new material and accretion of
smaller draplets at the island perimeter. In this early
growth stage, when two islands touch, they quickly coalesce
in to a larger island, and in doing so, wipe clean part of
the substrate which was covered. This wiping action
constantly creates gaps comparable to the island radius, and
results in ever-larger islands separated by gaps filled with
smaller islands nucleated since the last wipe.

At a later stage of growth, full coalescence no longer
occurs. Instead a partial coalescence of touching islands
occurs, resulting in wormy structures and little wiping.

To accurately mimic the film growth we have described ,
a model must have the following key features: an early stage
of growth and coalescence of compact islands, a crossover to
wormy structures and, under certain conditions, a very high
pc. The morphology of discontinuous films has been analyzed
in terms of percolation models. Various geometrical
exponents (fractal dimensionality, correlation length etc.)
have been measured in real films and are found to coincide
with values found in ZD pure percolation patterns. It may
therefore be natural to assume that simple lattice model
with short range attractive interactions will accurately
reproduce the morphology of real discontinuous films.

However, we have found that a large class of lattice models

49

fail to reproduce the essential phenomena of island
coalescence, worm growth, and high pc that occurs in the
growth of many thin metal films on room temperature
substrates.

The long and intertwined worms seen in Eignzg_1‘1 (lower
right) appear to be repelling each other. We will show that
a simple continuum model which incorporates the essential
feature of interrupted coalescence of large islands.and
which contains only short-ranged attractive interactions
accurately reproduces this apparent repulsion and high pc
which results from it.

Direct observation of growing metal films on warm
substrates has verified that full coalescence occurs when
the metal islands are small enough and partial coalescence
occurs when the islands reach a radius greater than a
critical value Re. In the regime that we are considering,
the substrate temperature is well below the droplet melting
temperature and the mechanism of coalescence is surface
diffusion (not bulk diffusion). The mechanism for the
interruption of full coalescence beyond a critical island
size is discussed in section 3.3. Here we consider Rc to be
determined by experiment.

The full coalescence stage has been previously studied

by Meakin and Family(Family and Meakin, 1988), and here we

50

extend their model to include the island-to-worm cross-over
and eventual percolation.

The MFM is defined as follows. All islands are assumed
to be circles (D-z model) or spherical caps (D-3 model).
Droplets of radius R0 rain down ballistically on a
substrate. If they land on a pre-existing droplet, a new
droplet of size

RD=R,°+R,° (3.1)
is generated. Here R1 is the radius of the pre-existing
droplet on the surface, while R2=Ro is the radius of the
droplet that is ballistically deposited onto the substrate.
D is the dimension of the depositing droplets (0:3 for
spherical caps, 082 for depositing circles). The center of
the new droplet is located at the center of mass of the old
droplet and the added droplet. If the new droplet now
overlaps an adjacent droplet on the surface, these two
droplets also coalesce according the rule (3.1). Complete
relaxation of the pre-existing droplets on the surface
occurs before the next new droplet is ballistically added to
the surface. According to this rule of coalescence, computer
simulation is easily carried out. Let the 2d coordinate
(x,y) stand for the position of a droplet which has the
shape of a circle. The original droplets have radius R0 and
are deposited randomly over system. The non-trivial physical

quantities are the coverage of film on the substrate and the

51

droplet size distribution. In the case 082, the total area
covered by droplets is independent of the size distribution,
therefore the coverage will be proportional to the number of
deposited droplets. For 083, the situation is different, the
coalescence rule conserves volume, the height of a spherical
cap is proportional to its radius. £1gn:§_1‘z shows the
droplet configurations of 30 spherical caps(Figure 3.2a) and
2D circles(Figure 3.2b) on square substrates. Indeed, they
have distinct size distributions which are statistically
shown in £1gg;g_;&1. £1ggzg_1‘1g shows the droplet size
distribution for the case D83. The distribution has two
maxima located near RsRo and RsRm, where RE is a length
scale that grows with total deposited droplet number. The
asymptotically bi-modal feature is not shared by D=2
droplets which have only one maximum near R=Ro(£iggzg_;‘1p).

The interrupted coalescence model (ICM) we introduce is
a natural generalization of the Meakin/Family Model(MFM) to
include a cutoff radius Rc above which overlapping islands
no longer coalesce. The ICM rule are thus:

a) If R1<Rc or R2<Rc droplets coalesce according to
equation (3.1). The new droplet is at center of mass of the
old droplets.

b) If R1>Rc and R2>Rc the droplets do not coalesce.

Rc sets the length scale at which the cross-over from the

droplet to wormy stage occurs. In the ICM we must also

52

 

":0 '0

5". .a‘.

' ‘ 5.. I .
0.. r «3“. .
O'O:.;.'P’...’91.

.fi'fi." 9

.- 03353-656:
.:’:::~'--;-f- "" - £91.”;
0393‘! d'. ' 3
’93"? '

 

 

 

‘0
.’

.'A.r‘.. e
‘ fié? =53. ‘ O o
’ '-°.‘.’:0;'§o (on
s“.
we.

. . .. . ‘
0.- ! o-
. tug} :2 k ‘ . ‘ 0 "0.9““
2.3. : ,_‘..e.-. ‘4': '
I . O" . o .. 3".)6. . o '0-

 

Figure 3.2: Droplet configurations on a square substrate using
MFMNmkin 8:. Family Model): a,b: D=3. c.d: D=2.

1030'.)

log(N,)(arb. unit)

53

 

lllllllllllllllllll

 

 

 

 

-i-

Ls

IlllllllllUlLllJlLll

 

 

 

1 l l l l l l I I l I l I 1 1

O
p

0
log(s)(arb. unit)

Figure 3.3: Droplets size distribution (a): 083, showing two
peaks,oneislocatedan=Ro;theotherincreaseswiththetotal
number of deposited primary droplets. (b): 082: the distribution
is monotonic, small size droplets have a large population.

54

define the radius R0 of the small droplets being
ballistically deposited onto the surface, and the edge
length L of the system. The important independent variables
in our simulations are thus the ratios Rc/Ro and L/Ro. We
consider the D=2 and D=3 cases, but restrict our attention
to 2-dimensional substrates (d=2), as the d=1 case has a
trivial PC. Finally, note that when Rc/Ro<1.0, no
coalescence occurs, and both the D=2 and D=3 models return
to the D=2 inverse Swiss cheese model (Pc=0.676i0.002).

Eiggzg 3.5 presents the growth morphology generated by
the ICM for the case Rc/R088, L/Ro-ZOO and D=3. It is seen
that in the early stages of growth, the model is very
similar to the MPH. As deposition continues, the large
droplets begin to overlap, until finally a percolation path
of the large droplets occurs. The configuration at
percolation (Eiggrg_ggfip) shows the high percolation
coverage generated by the ICM (in this configuration about
0.82).

The reason for the occurrence of high Pc is that the
initial droplet state is correlated, and that percolation
from this correlated state requires a higher coverage than
would random percolation of discs (inverse Swiss cheese).
The small droplets (R<Rc) do not lie on the final
percolation path, so the large droplets play the dominant

role. Given this, we can think of the final percolation path

55

1919130129.)“. 21:4.

 

 

 

Figure 3.4: [CM simulation with URo=200 and W880. (a):
‘Iheearlystage.Whenthelargestdropletsizeissmallerthanthe
criticalradiungthelCMsimuationreeultisidenticaltothe
MFM (b): Percolation stage: as more primary droplets deposit,
thosedropletswhoseradiiarelargerthankcfonnapercolation
structure.

$6

 

0.55 ‘ . .se . l - . . . T s . . e] -
0.50 _— i i f f
0.75 - i

0.70 -

L L A l l l l A | I I L L I

 

 

 

 

0.85 ‘
1 2 3 4
Rc/Ro

Figure 3.5: Percolation coverage as a function of Rc/Ro for D=3:
The percolation coverage is calculated fiorn averaging 100
configurationsforcachpoint'l‘heerrorbarsstandforthe
standard deviation. As Reiko increases from 1, Pc has a bigjump,
which contributes to the high coverage.

57

as evolving from a set of repelling discs (the large
droplets) with packing fraction about 0.57 (as measured on
configurations such as Eigg;g_1‘5§). If we now increase the
radius of these random discs (at a rate proportional to
their original radius - to model the cross-section available
to the ballistically depositing droplets), we naturally find
a high PC. If the original droplets lie on a triangular
lattice Pc=0.907, while if the initial state is like random
packing, we expect 0.82<Pc<0.89. The point is that starting
from a distribution of repelling droplets always leads to
high PC, and given an initial droplet configuration such as
Eiggzg_;&5g, it is not surprising that the ICM model does
too.

It takes the MFM many iterations to settle into a self
similar evolution. Prior to this equilibration time, the
droplets are not as strongly correlated. We would thus
expect that small Rc/Ro implies lower PC. This expectation
is supported by the data of Eiggzg_1&§, which shows results
for Pc as a function of Rc/Ro for D=3.

We have also studied the D=2 ICM model. Although the
MPH in D=2 leads to a very different distribution of droplet
sizes (bi-model in 3-D, monotonically decreasing in 2-D),
the final percolation coverage in D=2 is still enhanced in

the ICM. This is illustrated in Eiggrg 3.§ which shows two

 

 

 

 

( o.

(b)

 

 

Figure 3.6: ICM simulation for D=2 with UR0=200 and
Rc/Ro=8.0. (a) early stage: similarly the early stage is identical to
the MFM case. Droplets are isolated. (b) percolation stage: the
introduction of the critical radius naturally leads to the
percolation structure.

0.85

0.80

0.75

0.70

0.65

59

 

 

 

chRo

Figure 3.7: Percolation coverage as a function of lie/Ro for D=2 :
Here the percolation coverage is calculated from averaging 200
samplesatpercolationforeachpoint. ltisalsoscenthatasRcho
increases from 1, Pc has a bigjump, but smallerthanthat in D=3.

 

60

examples of the D=2 morphology and in Eiggzg_;&l which
presents Pc(Rc/Ro)

An important feature of the dependence of Pc on Rc/Ro,
in both D=3 and D=2 (see Eiguzg§_1gi_gng_1gl) is that there
is a first order jump in P6 as Rc/Ro—n. At Rc/Ro=1.0-,
Pc=0.676i0.002 in both D=2 and D=3, while we find that
Pc(1.0+)=0.745i0.005 (D=3) and Pc(l.0+)=0.740i0.005 (D=2).
These non-analytic changes in Pc at Rc/RO=I.O are physically
reasonable, because, as soon as Rc/Ro>1, all the droplets
that are being ballistically deposited onto the surface will
coalesce with any pre-existing droplet with which they
overlap i.e. an infinite sequence of coalescence events
occur in an infinite system. When Rc/Ro=l.0- however, no
coalescence at all occurs.

Since there is a jump singularity in Pc at Rc/Ro=1.0,
we must consider whether the universality class of the
percolation geometry is altered for all Rc/Ro>1. We have
tested the universality class of the percolation problem at
Pc in the ICU through calculating the cluster size
distribution, and find that, in the medium size range, the
distribution complies with the scaling law very well. figure
lgfi shows the cluster size distribution at the percolation
point. According to the scaling theory of percolation, the
size distribution na near the percolation point has the

scaling form:

61

 

    
  

 

 

 

 

0° Tfi I I I i T T ITII I"! Till "II "II" I I I I
1 °°5 I I I I I l l "IT I n.-
l. I
5000 -- (a):D-3 -
2000 - ' o s —
z' 1000 :- —.
500 '— —'
- cl
200 '— —‘
101ml= ~
1: 10 20 E
5000 :- (b):D-=2 -:
c: 1000 .— -:
a 500 :— -‘
. l. 2
.0
h I- "l
3
50 :—
C
10 l l l l l l l l l
1 2 5 10 20 50 100
s(arb. unit)

Figure 3.8: Cluster size distribution from KM (11) D=3, in the
modiurnsizerangeJhelog-logplotisnarlyastraightlincboth
the small and large cluster size fall ofi‘the line significantly. When
sissmall,thescalingbehavi0rfails;whcnsistoolargethesize
efi‘cct has major influence on the scaling behavior.

62

".(p)=3"f[(P-P.)S"] (3.2)
where t and o are critical exponents. The precise form of
the scaling function f(x) has be to determined by computer
experiments and other numerical methods. Other critical
exponents are to be expressed using I and 0. Here we omit

the derivation and give the results only(Stauffer,1985):

a=2-(r-l)/a (3.3a)
fl=(r-2)/a (3.3b)
y=(3—z)/a (3.3c)

t is measured from the slope of a log-log plot of n,, y can
be determined from the measurement of average cluster size.
The values are listed in Igblg_3‘1 which shows the

comparison with conventional percolation results(Stauffer,

 

 

 

 

 

 

1985).

Table 3.1

2d 20 on 2d 30 on 2d
a -2/3 -0.71i0.08 -0.70i0.08
B 5/36 0.18:0.08 0.1510.08
1 43/18 2.35:0.04 2.40:0.04
t 187/91 2.07:0.04 2.06:0.04
0 36/91 0.39:0.04 0.39i0.04

 

 

 

 

 

 

It is seen that within a comfortable tolerance(less than
10%), the critical exponents in the ICM model are matched

with the results from percolation theory. Therefore, the ICM

63

percolation process appears to remain in the same
universality class as the uncorrelated case.

We note that the wormy structures of Eiguze 3.1 do not
at first sight resemble the percolation structure of Figure
3.4b. The main reason for this is that sharp grooves occur
at the intersection of the circles in £1ggrg_1‘1p, but do
not occur in Eigg:§_;‘1. In a real film these grooves would
round out to a radius of curvature RC, producing the wormy
appearance so characteristic of these films. Unfortunately,
we have found no efficient method for simulating this
rounding process in random continua shown in Eigure 3.52.

In summary, we have introduced a simple continuum model
of film growth. Our model is an extension of one by Meakin
and Family and includes a cross-over from an early island
stage to a later wormlike stage. This model provides a
natural explanation for the high Pc observed in many thin

metal films grown on insulating substrates.

 

As stated in last section, in many film/substrate

systems, the depositing material does not wet the substrate,
so that in the initial stages of growth, distinct islands of
the deposited material form (Volmer-Weber growth)(Pashley

etc., 1964). As deposition continues, these distinct islands

reach a critical size Rc(Rc-dependent on substrate

64

temperature T, deposition rate and anneal time amongst other
things) at which they connect together to form a percolating
structure, and further deposition soon leads to a continuous
thin film. The case of layer by layer growth is the
exception rather than the rule, as films may be rough either
due to non-wetting, strain mismatch, due to kinetic
roughening, or due to a combination of these processes. The
early stages of film growth, and in particular the value of
Re are important in controlling the final morphology of
polycrystalline films.

A detailed experimental study of the Re dependence on
deposition rate t and substrate temperature T for Pb has
been performed on very smooth and clean amorphous Sioz
substrates(0ubson and Jeffers 1994). To explain this data, a
kinetic freezing model has been developed and provides a
semi-quantitative explanation of the deposition rate and
substrate temperature dependence of Re. In the kinetic
freezing model, it is considered that the complete
coalescence time for two spherical droplets is proportional
to the fourth power of radius, and thus the process of
complete coalescence is overtaken by deposition. However, we
also noted that the kinetic freezing model predicts a
stronger than observed dependence of Rc on deposition rate
(Jeffers etc., 1994), and suggested that this reduction in

coalescence may be due to substrate inhomogeneity pinning

65

the islands in meta-stable states. Suppose that the surface
energy is isotropic, the gain in surface energy in forming a
compact shape from an elongated single crystal island scales
as R2, while it will be seen that the pinning energy scales
as R3/3. It would thus seem unlikely that pinning can
produce metastability of elongated single crystal metal
islands. In contrast, an elongated metal island containing
one or more grain boundaries in its interior is easily
pinned in a meta-stable state. The combination of grain
boundaries and substrate inhomogeneity thus provides an
important mechanism for strong substrate pinning effects on
the growth morphology of polycrystalline films on non-
wetting substrates.

Eiggzg_1‘2, gives experimental results for TC(R), the
critical temperature below which no complete coalescence
takes place at fixed deposition rate. It is obvious that
Tc(R) increases with R. The experimental data for the
Pb/Si02 system is shown as the solid dots in Eiggrg 3.2.
Islands of size R>RC(T), touch and form a contact, but do
not coalesce, and this region has been labeled sintering as
it is similar to neck formation during the initial stage of
sintering of granular aggregates. The dotted line is a
schematic of experimental results on the bulk melting of
lead, which has a size dependence TM(R)=T,,,(00)(1—c/R), where c

is related to the densities and surface energies of the

suasrnArE TEMPERATURE (K)

66

 

 

 

     
 

 

 

65° 1 1 I I I I r
BULK mourn suut uEero: 001 it

”I” """"""""""""""""""""""""""""""""""" “

550 L euut soup —
FULL COALESCENCE

- —a_-
450 .. -

* SINTERING ‘

(NECK FILLING)
350 - -
250 4 l 1 l 1 J 1
0.0 0.5 1 .0 1.5 2.0
CRITICAL RADIUS Re (um)

Figure 3. 9. Experimental measurements for T¢(R) The upper
solidlineisthesizedepcndenceofmeltingtemperatures; the
dottcdlineisthecurvefromexperimcntaldata;theareabenveen
dianisdiecoalescawercgiominwhichsurfacedifi‘usionisthe
dominantfactorforcoalesccnce.

67

liquid and solid phases, and to the latent heat of the
liquid-solid transition. Typically metal clusters have to be
of radius of order Snm before a significant (of order 10%)
depression of bulk melting occurs(Buffat,1976). Between the
bulk melting curve and the curve Tc(R), there is a broad
regime in which solid metal islands coalesce, predominantly
by surface diffusion. In this regime, the island size
distribution can be described by breath figure models. As
the island size increases during growth, the film crosses
the curve RC(T). At this point the metal islands do not
fully coalesce into compact shapes, elongated metal islands
are produced, and the film reaches the percolation point
shortly afterwards.

Since the substrates we are considering are amorphous,
separated metal islands grow with different crystallographic
orientations. Thus when they first overlap, there is a grain
boundary at their intersection. This grain boundary is
eliminated during the coalescence process, as evidenced by
the fact that for droplets R<<RC(T), each grain is a single
crystal. However, when elongated metal islands first form,
they usually contain a grain boundary in their interior,
suggesting that grain boundaries play an important role in
the energetics of coalescence.

When a grain boundary intersects the free surface of a

solid, it forms a grain boundary groove. The equilibrium

68

 

 

 

 

 

 

 

(a)
\\\\W \
“““ieswix ““W ‘1“““ \\
I)! W
(b) . .
a l3
e > i <
R1 R2

Figure 3.10: Two touched droplets with an interior grain
bnoundary.’l'helargerdr0plethasradiusR¢andthesrnallerhas
radiusR,.(a)sh0wingthedefinitonofv(b)showingthe
definitionofcandfl.

69

grain boundary groove angle, w, is related to the surface
energy of the metal, 7, and the grain boundary energy, yb.
In the case of isotropic surface energies, the geometry is
as depicted in £1ggzg_1glga, which leads to the expression,
71. =27008¢ (3.4)

' For the purposes of illustration, we specialize to the case
of the wetting angle 9-90°, and with the variables as
defined in Eigg;g_;‘1g, the sum of the surface and grain
boundary energies is given by:

E = n7[R,2(1 +cosa)+R22(1+cosfl)]+ ngn—cosm- w] (3 .5)
where

Rb=Rlsin a/sin(a- V!) (3.5a)

The angles a and B are as defined in zigg;§_gglgb. Here, we
have used the fact that the grain boundary is a segment of
the surface of a sphere, which bisects the groove at the
triple point. It is assumed that the grain boundary energy
is isotropic. Although the surface energy is certainly
anisotropic, as evidenced by the anisotropic shape and
significant flattening of parts of equilibrium crystals,
this should not change our prediction for the scaling
behavior of Tc(R). In fact the only thing that is important
for our argument to hold is that the presence of the grain
boundary lead to unstable equilibrium of an elongated metal

island. This feature should also be true for appropriately

7O

chosen anisotropic equilibrium and non-equilibrium crystals
shapes.

To ensure that the volume, V, of metal is conserved
during coalescence, we impose the condition,

V=(7r/6)[R,’(2+3cosa—cos’a)+R§(2+3cosfl-cos3fl)]=const. (3.6)

For 6<(90-w), the grain boundary groove tends to grow all
the way down to the substrate(Hillins,1957), and hence
separate the elongated island into two single crystal
grains. If 9>(90-w), there is a finite groove depth, and
the energy by equation(3.2) is a monotonically decreasing
function of the size of the larger metal island R2(see
£1gg:g_;&;1). In particular, when the two metal islands are
nearly the same size R/Rzzl, there is a very small slope
in the surface energy curve and hence the driving force
toward coalescence is also very small. It ought to be
pointed out that the small slope is due to the presence of
the grain boundary. In the absence of a grain boundary,
there is a steep slope in the analogous surface energy plot
of elongated metal islands. Since the energy curve(£igure
1‘11) is flat for nearly equi-sized islands, relatively weak
contact line pinning, by for example surface irregularities
or impurities, is able to prevent the metal islands from
coalescing. Once a pair of islands is pinned in a meta-
stable state, islands which are of similar or larger size,

which attach to the pinned pair further stabilize the

total surface energy Ely

71

 

8.1 -
8.0 b
7.8 ’
7.8 b
7.7 .
7.. l

7.5-

 

 

 

 

 

 

 

 

V=70° m _

1 1 1 '1

 

 

J

 

1.00

I 1 1
1 .02 1.04 1 .06 1 .08 1 .1 0
radius 01 larger Island R,

Figure 3.11: Surfice energy as a function ofthe radius ofthe
larger’dropletR1.WehavetakenR,-1.Whenlt,iscloseto l,the
surfaceeneryisveryflatwhichmeansasmallpinningforcecan
stopthecoalescenceprocess.

72

structure. In this way, grains aggregate to form the
percolation structure characteristic of the morphology of
films grown on non-wetting substrates. Although the early
coalescence stage does have a strong effect on film
morphology, the large scale geometry at percolation is
insensitive to this short range correlation, as expected
from universality arguments.

We now concentrate on the case 9>(90-w) which is
applicable to many metal/insulating substrate systems, where
the wetting angle is typically quite large, but 90-w~10o for
high angle grain boundaries (where nfq/3). We develop a
scaling argument to estimate the temperature at which
thermal fluctuations can activate a pinned elongated metal
island from a meta-stable state. We build on the analysis of
contact angle hysteresis in fluid flow through capillaries
(Joanny and de Gennes,1984). Here we calculate the energy
due to the slight deformation n(x) of a contact line. The
surface profile of the liquid in the absence of bulk
external f:;ld‘;atisfies Laplace equation:

2 z

Est-5?: (3.7)

with the boundary condition that z vanishes on the contact
line: z(x, r)(x)) = O . For 71(x) = O,z(x,y) = ysin 6, which corresponds to
the ideal case. The solution of equation (3.7) can be

expressed as:

73

 

. sinfim , ,
z(x,y)=ysrn G-Tirflx)y2+(:-x.)2dx (3.3)

The incremental energy due to the presence of n is expressed

as

E. =-;- I drdyrtIVzY—(sin 6021 (3.9)

where y is the surface tension of liquid. Inserting equation

(3.8) into (3.9), and integrating over y, one has

. d
E.= 78m BIdr'fU) (3.10)
0

 

d

Here 1/d is a low cutoff wavenumber. Defining

. 1"
2r zgjrflxyir (3. 11)
0

as the deformation magnitude, the energy increment becomes:

Ed=%klf (3.12)

Where k==ymn26,independent of system size. Any deformation
of the contact line will cost energy. In reality, the
substrate is far away from perfect smoothness and
homogeneity, therefore a deformation may be energetically
favorable. Suppose the correlation length of random

potential is smaller than the deformation size n, the energy

gain is proportional to square root of deformed area JnL(

note that it is.JnL not nL due to randomness)(Robbins and

Joanny,1987) where L is the contact line length
Eh=-h./r)l. (3.13)

74

where h is constant depending on pinning strength and liquid
surface tension. Thus the total energy change due to the

slight deformation of contact line is:

1
E(U.L)=5kfl’-NUL (3.14)
The minimization of equation (3.14) gives the deformation
magnitude
Lhz 1/3
—- 3.15
"4“[41‘2] ( )

and the corresponding energy gain:

10
5,0.) = 4.75116?) L”3 (3.16)

while L=2nR, it thus follows that the pinning energy scales
as R3/3, as mentioned previously.

Although this argument was developed for liquids, it
applies for solids provided a significant fraction of the
exposed faces of the crystal are above their roughening
temperature. If the temperature is sufficiently high,
thermal fluctuations depin the contact line. We estimate the
depinning transition by comparing the thermally induced
contact line wandering with the disorder induced wall
wandering of equation (3.15). To calculate the typical
amplitude of thermal wall wandering, we consider the height-
height correlation expression

C(r) = ([z(r)z(0)]) (3 . 17)
The Hamiltonian related to the surface profile is roughly

expressed as(Week,1980):

75

m=fidxayBy(V-:)2] (3.18)

where y is surface tension. From the Boltzmann distribution,
the correlation function becomes:

G(r) = ([z(r)z(0)]) = E i”— (3 . 19)

2
‘7

Turning the sum over q into an integral and introducing a
unit lower cutoff scale(lattice spacing), the thermal
fluctuations leads to the wandering amplitude

17, z(kBT/Znyln R)"2 /sin 0 (3.20)
Equating nT and flat

k,7;(R)ch2’3/lnR (3.21)
with c given by,

cz27rysin «trill/21:2)”3 (3.22)
Notice that k in equation (3.21) relates to y and 9 through
k=2)'sin 9, thus the coefficient c is inversely proportional
to (sin9)'2/3. The smaller 9 is, the more effective is the
pinning, and hence the larger Tc(R) is (of course this
breaks down if 6>(90-w)is violated). As expected, larger h
implies larger c, which leads to larger TC(R). These general
conclusions hold even when more complex forms of thermal
Hamiltonian are used, although the R dependence in equation
(3.7) is clearly altered. Although based on several
oversimplified assumptions, equations (3.20)-(3.22) provide
a useful guide to the expected effect of substrate

inhomogeneity on growth morphology.

76

So far, the experimental results are for Pb on very
smooth and clean amorphous 3102, and even there, there are
clear indications that both pinning and kinetic freezing
effects are important in determining film morphology. It
would be useful, although difficult, to develop a theory
addressing the combined effect of kinetic freezing and
substrate pinning. Another useful approach is to try to
experimentally isolate these two effects by varying
substrate roughness, substrate temperature and deposition
rate. For sufficiently strong substrate disorder, it is
probable that a strong pinning theory would be needed

instead of the weak pinning result used here.

CHAPTER FOUR
MOLECULAR DYNAMICS STUDY OF METAL CLUSTERS

W

In this section, I will review some fundamentals of
statistical mechanics and thermodynamics that are relevant
to our analysis of computer simulations. We assume all
particles are identical, have three spatial degrees of
freedom, and obey classical mechanics defined by a total
Hamiltonian flfl In the canonical ensemble, the number of
particles N, the volume V and the temperature are held
fixed. The motion state of the system is represented with
generalized coordinates r=(r1,...rN) and conjugate momenta
p=(p1,...pN). Generally, fl’is a function of r and p. The

probability density at (r,p) in phase space is of the form
‘P(N.V. T) = CXPl-WUJH (4 - 1a)

 

l
Nlh3”Z(N,V, T)

with
l
AHh“

 

2(N.V.T)= fidvdpepr-Mrmn (4.1»)

The fundamental relationship between statistical mechanics
and thermodynamics is given by

T(N,V,T)=—lenZ(N,V,T) (4.2)
where T is Helmoltz free energy and is minimized at
equilibrium. The Hamiltonian fl'is a sum of kinetic energy

71p) and potential energy Ikr) for most systems we are

78

discussing: ififl4tt Therefore, the momentum part of the

partition function Z can be integrated out:

 

__ Q(N,V,T)
Z(N,V,T)-————N!A3N (4.3)
where Q and A are given by:
Q = j epr-flWerr
4.3
(2ka )1/2 ( a)
A== 2
h
The canonical average of a variable is given by:
(A) = j Aexp(—fi}{)/Z(N,V, T) (4.4)

The average of T'thus is easily obtained: CT)=3A%T/2.
Molecular Dynamics simulation calculates the positions and
velocities of particles in the systems from the Hamiltonian
described above. The equations of motion of each particle

comply with Newtonian Dynamics:

gEv--——- /nr (4 5a)
dt asp. 9" '

. ax av
5%“?“7 “-5”

In traditional molecular dynamics, the total energy E for a
fixed number of atoms in a fixed volume V is conserved as
the dynamics of system evolves in time, and the time average
of any property is an approximate measurement of the
microcanonical ensemble average of that property for a
thermodynamic state of N,V,E. However, for certain
applications, it may be desirable to perform dynamic

simulations at constant temperature. Several artificial

79

inventions are used to achieve this goal and we will focus
on the Extended System Molecular Dynamics(ESMD) approach
(Nose,1984; Jellinek and Berry,1989; Anderson, H.C.,1980)
which allows the system to thermally contact with a heat
reservoir represented by an additional degree of freedom s.
To realize the objective that the average kinetic energy is

fixed at equilibrium, the Hamiltonian is postulated to be:

9f=2pfl2msz+'V(r,)+p,2/2Q+nglns (4.6)

P. is the conjugate momentum of s; Q is a parameter of the
heat reservoir and behaves as a mass for the motion of s; kT
has its regular meaning; it will be shown that g is
essentially equal to the degree of freedom of the physical
system. It is also assumed that the Hamiltonian equations of

motion of r and s are valid:

 

Enhflbfl

dt 51-,- a

dr,_o‘t1-[_ pi

dt - d),- -m,sz (4.7a)

and

ch 65 s t

(4.710)
é:__£{_=m
d1 40.. Q

It is easy to prove that 5{ is conserved through combining
equations (4.6) and (4.7). Therefore, the micro-ensemble

partition function can be expressed as:

 

the

10

vi

80

z=jdp,jdsjdpjdrx5{a{-E] (4.8)
where E is total energy. If we use p'=p/s, r'=r, and set
g-3N+1, the degree of freedom of the system, 2 yields the
following expression after integration over p, and s:

z=cfidp'dr'exp[-fi(q'+a/)] (4.9)
where C is a constant. This is exactly the canonical
ensemble partition function. The MD simulation is to solve
the equations of motion (4.7) for appropriate initial

conditions.

- . -

One very important aspect of molecular dynamics
simulation is to have a good potential describing the
interactions of atoms in the system. The Lennard-Jones
potential is famous and successful when applied to rare
gases, e.g. Ar and Xe. But it fails to predict and explain
many experimental data of metals, especially the noble and
transition metals. For instance, the Cauchy relation of the
elastic constants C12/C4¢=1(Landau,1984) is valid for any
pairwise central potential, of course including the Lennard-
Jones(LJ) potential while the experimental data for this
ratio is much higher: 3.7 for gold and 1.5-3.3 for
transition metals(Ercolessi etc.,1988). Another failure of
the LJ potential and other two-body interactions is that

they predict the cohesive energy BC, the energy needed to

81

remove an atom from the bulk lattice is essentially equal to
the vacancy formation energy Ev; while the experimental data
show that EV/Ec80.25 for gold and about 0.35 for transition
metals. In a two-body system, the melting temperature is
usually near 0.11%; in gold, on the other hand T_~O.O3Ec/k8.
Some of those problems may be cured by adding new terms
which depend on the total volume of the system in the
potential. Therefore, it is clear that, to describe metals,
an adequate model must include many body contributions.

In metals, the conducting electrons are loose and their
collective Coulomb interaction with ions is the force which
maitains the lattice structure. For transition and noble
metals, the conduction electrons are d-electrons and d-
electrons. Generally the s-electrons form a very broad
energy band which overlaps with the narrower s-electron
band. The following calculation based upon the work by
J.Friedel(Friedel, 1969) and Gupta(Gupta, 1981) is to sketch
an interaction form starting from the tight binding model.
First, one considers d-electrons which have 10 atomic
states. For the sake of simplicity, the density of states is
represented by a rectangular function of width W and height
10/W. The Fermi energy E, for a metal with number 2 of d-
electrons is related to the bandwidth through:

EF=W(Z-S)/10 (4.10)

82

The zero energy is at the center of the band. The cohesive

energy Ed from the d electrons is given by:

E
Ed: [545% (4.11)

47/2
where p(E), the density of state has been assumed to be

constant. The integral then gives:

W’
E =—Z lO-Z 4.12
. 20 ( ) ( 1

In the tight-binding model, the bandwidth W is determined by

the overlap integrals. The second moment of p(E):

m2
1

l
=_ 152 E)dE=-W2 4.13
#2 1041/2 A 6 ( )

On the other hand, 0; can be expressed as:

u. =Z(kI(H-E.)zlk) (4.14)

h
where H is Hamiltonian of the d-electrons and BC, the atomic
energy level at a selected site i can be taken to zero. k
labels the Block wavevectors. H2 is expressed in real space,
reading:
14. = (ilH’li)= Zl<iIHlj)l’ (4.15)
2
where |i> denotes an atomic wavefunction of a d-electron at
site i. If only nearest neighbor overlap is taken into

account,

m=ZfW> MAM

131'

Here B(rij) is the hopping integral and ri is the distance

i

83

between sites i and j. Combining all these equations and

eliminating W, one obtains:

Ed: 232(5)] (4.17)

131

with C=(\/6/20)Z(10-Z). The hopping integral can be taken to

be of exponential form:

Ar) = 16.91"" (4 . 1a)
where a is the bulk interatomic separation at equilibrium.
The equation above gives the d-band contribution to the
cohesive energy. The stability of the lattice requires that
there be a countervailing short range repulsive force, which
is provided by the compression of the free electron gas, in
the model, the s electrons. As the atoms are brought
together to form a solid, the free electron contribution to
the cohesive energy is at first attractive since the valence
charge density is expelled into a region of more attractive
potential. However as the interatomic separation is further
decreased, the potential reaches its minimum and the
repulsive force begins to dominate. According to this

picture, the short range repulsive force can be written as:

E, = 82 e’M‘“ (4 . 19)

jar
and the parameter p characterizes the interaction range. The
total energy Ec=Ed+E8 ought to be minimized at r=a.

Consequently, the following expression can be obtained:

84

-1 - - 1/2- q - —
Isl-2V. [Zexpl 240.,- a)]] WZexpl p(r., a)] (4.20)

131 121

where z is the number of nearest neighbors. For a specific
metal, the parameters p, q, and V0 can be determined by
experiments. While doing computer simulations, one just
needs the reduced expressions and any quantity in the final
results can be converted into real units by dimensional
analysis. To do that, we set a=1, Vb=1, thus the total

energy of a metal system with N atoms can be written as:

”2
N
V = —%Z [Ze‘mfl’] - AZe'M‘” (4 . 21)

F. p.
where A=ql pf . This energy form cannot be modeled by
pairwise forces because a two-body scheme implies a linear
dependence of energy of an atom upon its coordination;
Because of the presence of the square root term in the
expression, the energy is nonlinearly dependent on its
coordination. The potential has been used by many authors
(Bulgac,1992; Garson and Jellinek,1992) and proved to be
able to describe the bulk and surface properties of a
variety of metals well. Next, I will apply the potential to
gold clusters and study their equilibrium and dynamic

properties.

85

. - d
The notion that small clusters might exhibit

distinguishable solid-like and liquid-like phases arose from
classical isoenergetic molecular dynamics(MD) simulations of
the 1970's(Briant and Burton,1975), largely of Argon
represented by pairwise Lennard-Jones interaction. Results
'of simulations indicated clearly that solid-like behavior
could be anticipated at sufficiently low energies and that
liquid-like behavior would appear at high energies( Berry,
1987; Jellinek etc.,1986). The typical curve of the
temperature versus energy is shown in Figure 4.1. This type
of behavior appeared for clusters with N>6 and was
interpreted as indicative of a first-order phase transition.
More recently, different potentials describing various
systems such as metal and carbon clusters have produced
similar results. The striking evidences of those simulations
inevitably raised the question of what the physical basis is
for the existence of solid and liquid phases in finite
systems, especially small clusters. The following argument
based on the thermodynamic point of view shows that it is
possible for small clusters to have a fairly sharp solid-
liquid phase transition. A solid-like phase can be
intuitively identified with a state of rigidity and liquid-
like phase with a state of floppiness. Hypothetically,

suppose an exact Hamiltonian describing a cluster could be

86

 

-0.80 1111111111111711111
l l l l

-0.85 — __
-0.90 —' —
-0.95 — -

 

 

l l l l I l l l 1 l l l l 1 l l l l I l l l l I

.00
0.0000 0.0025 0.0050 0.0075 0.0100 0.0125
'1'

 

Figure 4.1: A solidoliquid phase transition schermtic plot: total
energy v.s. temperature: The calculation was done for cluster
N-13 using both constant temperature and constant energy
moleculardynanucs.1‘hetw0niethodyieldsthesanieresult1'he
solid-hquidphasetransitionissignificdbydiepresencecfanvist
rcgrononthecurve.

87

solved, one would have the exact statistical properties of
the cluster and then one would know the cause of the phase
transition. In a theoretical study, most clusters are
customarily described by effective Hamiltonians that are
approximate in two ways. They are based on the adiabatic or
Born-Oppenheimer approximation which supposes separability
of the slow nuclear from fast electronic motion, and on the
assumption that the amplitudes of vibration are small
compared with the interatomic distances, the frequencies are
high compared with rotational frequencies and the moments of
initia are therefore nearly constant. The energy levels of a
cluster with the above approximation are described as the
electronic excitations plus vibrational spectra which are
surrounded by more detailed rotational energies. The density
of states D(E) of a system in such a rigid phase is higher
than in a floppy phase in the lowest range of energies. The
densities of states of both the solid-like and liquid-like
limits increase with energy, but D(E) of the liquid
increases much faster so that at sufficiently high energies,
the density of states of the liquid-like form becomes much
larger than that of the solid-like form. Here we introduce
an order parameter-like quantity 0 to represent the
floppiness of a cluster, when ¢=1 the cluster is in the
completely floppy state while ¢=0 is the limit for the

solid-like phase. Consequently, at low energies, the

88

partition function Z(T,¢) mainly has a contribution from low
energy states and thus is a decreasing function of 0. On the
other hand, at high energies, Z(T,¢) is an increasing
function of 0. The free energy F(T,¢) has only one minimum
at 0:0 at low temperatures and one minimum ¢=1 at high
temperatures. As the temperature is raised from the low
energy where only the solid is stable, D(E) near 0:1
increases faster than that near ¢=0, so that Z(T,0)
increases slower with T than Z(T,1), and then F(T,0)
decreases slower with T than F(T,1). As T increases, it
reaches a value Tf at which the free energy F(T,¢) develops
a slope of zero at some value of 0 at or near 1; that is
d7(T,¢)/0"¢l,=,,=0. At temperatures above Tf, F(T,¢) has two
minima. As T increases still further, it reaches a value Teq
for which the free energies at the two minima are equal.
Because F(T,¢) becomes monotonically decreasing at
sufficiently high T, and is a smooth function of both T and
0, there must be a temperature Th at which the minimum near
¢=0 turns into a point of zero slope; above Th only the
liquid-like form is thermodynamically stable. Eigure 4.2
schematically shows how F(T,¢) changes with T and 0. So far,
we have used the properties of the densities of states to
argue about how the free energy of a system changes as
temperature increases. It should be pointed out that, in

general, the range Tf<T<Th is the transition region. In the

F(¢.T)

89

 

T<<Tm

 

 

T=Tf

 

Tori

 

 

T=Th

 

 

T>>Tm

 

 

¢

Figure 4.2: A schematic plot of free energy F as a function of
floppiness ¢ at changing temperature. At low temperatures, F has
onlyomminimumatoio,thesolidstate;atT=T,-,Fisflatato
=1;whenfi1rtherincreasingT,Fhavetw0minimaat¢=0and¢
=1; at T-Teq, the values at thwe two points are equal. which is
coexistingphase;atT-=Th,Fhasonlyoneminimumat¢=l, the
liquid state, but F is flat near ¢=O; At high temperatures, F has
onemrmmum.

90

following discussion, the transition temperature( melting
temperature) is defined from our MD calculation results.
Both experimental data and simulations show that the
melting temperature Tm is dependent on cluster size N. The
experimental data are summarized in £1gnre_5r3 for Gold
(Buffat and Borel,1976; Borel,1981). The dots represent the
experimental data and the solid curve is the fitting

function which is of the following empirical expression:

T, =T,(oo)(1-5r‘-) (4.22)

where the bulk melting temperature, Tmon)=1338K for gold.
The following argument presents a verification of the above
relationship from the thermodynamic point of view .
Consider a solid gold spherical particle of radius r.
Its total free energy includes the surface term and bulk
term, reading: F,=Mp+A7. M and A are total mass and
surface area respectively. 7 and u are surface tension and
chemical potential. Suppose the same amount of gold is in
the liquid state, the free energy is rewritten as:
F, =Mp'+A'7' Notice that the total mass remains unchanged
while the surface area is altered because of the density
difference of liquid and solid gold. At the melting
temperature Tm, F;=F,, which leads to M(p-;1')=Ar—A'r'.
Assume that Tm is not very far away from the bulk melting
temperature Tmow), u could be expanded and

p(T,,,)-,u'(T,,,)=p,-p,-AT(S,-S,) where p,=,u3 at Tm(oo) and S is the

91

mo. 0qu --.

1200 .

1000 ..

 

 

Figure 4.3: Experimental data of Tu(r): the solid dots and
trianglesdenotesexperirnentaldatawhilethesolidcurveisthe
semi-empenealfittothosedataandhastheanalytical form
expressed by equation (4.22).

92

entropy of unit mass. S,-S,=~L/7;0w) at the bulk melting,
where L is latent heat of unit mass. AT=Tm(oo)-T,,,. The
geometry of liquid and solid sphere has the following

relation:

. 2 _ 21 2,3
A7-A7'=47'(’.7 ’17): 3 [7-71[£’2) ] (4.23)

 

M Ami/3 m. p.

Combining all these equations, one has:
25
AT 3 ,‘p]
= 7‘7"' (4.24)
Tn(w) rpL[ [pl ]

Here we have dropped the subscript of solid. The equation is

 

consistent with the empirical result equation (4.22) with rc
given by r.=(3/pL)[7-r'(p/p.)”’].

To verify the occurrence of a phase transition utilizing
MD simulation, we studied the caloric curve i.e. the
relationship between the internal energy E and the
temperature T which is proportional to the average of the
kinetic energy. In the ESMD( extended system molecular
dynamics) simulation T is an input parameter. Any computer
simulation is a speed and accuracy trade-off. The minimum
requirement for MD simulation is to satisfy the conservation
laws(total momentum, total angular momentum and total
energy) in the predefined tolerance ranges. The computer
program is printed out in Appgngix_5 in which we utilized
the predictor-corrector algorithm to solve equation (4.7).
The first step of the MD simulation is to choose the atom-

atom interaction. The potential we have used in describing

93

gold clusters is equation (4.20) discussed in figgrign_1rz.
The parameters for gold are a-0.288nm, Vb=3.81eV, p=10.15
and q=4.13. In addition, the mass of the atom has been taken
to be one so that the force on an atom is equal to the
acceleration at that time. The conversion from a reduced

unit to real time unit is summarized in Table 4.1:

 

 

 

 

 

Table 4.1
conversion red. unit real unit
time t 1.0 2. 12x10'13s
length a 1.0 0.288nm
energy E 1.0 3.81eV
temperature T 0.01 441K

 

 

 

 

 

The time step St has been chosen from 0.002 to 0.02,
depending on the temperature; high temperatures require
small time steps and low temperatures could have relatively
large fit for the same accuracy requirement. But in our
calculations 6t never exceeds 0.02 in reduced units, which
secures the deviation of total Hamiltonian H in
equation(4.6) within 10". The caloric curves are obtained
for four different cluster sizes: N=55, 177, 381, 725 which
roughly corresponds to spheres of radii r=2,3,4,5
respectively. The initial conditions have been set to be of
total zero momentum and angular momentum, preventing
clusters from drifting or rotating. Clusters are set in a

vacuum(free boundary) and low temperature environment at the

94

beginning. To observe the correct equilibrium quantities,
the system must reach equilibrium before any quantity is
measured. To achieve that, the system is run about 10,000
(0.1ns) time steps and a set of physical quantities(
typically internal energy E and kinetic energy K) are
attained. A subsequent 10,000 time steps are taken to
compare statistical results of those quantities with
previous ones. Typically 10,000-100,000(1ns) time steps are
entailed for the system to reach equilibrium. The total
internal energy E is calculated from averaging over another
100,000 time steps. The results are illustrated in Eiggrg
5‘1. Each curve has a sharp increment in E as temperature
increases past a certain value as indicated in the figure.
It is very noticeable that the increment in E and the
temperature indicated in Eiggrg 5.5 become greater as system
size increases. A more obvious definition of the transition
temperature is from the specific heat curves CV(T) that are
easily found from the caloric curves(Figure 4.4a). Each
curve of CV(T) has a maximum defined as the melting
temperature 7; which moves to the right as N, the system
size increases. Also the peak of Cv sharpens and approaches
to the bulk limit where Cv diverges at melting. Eiggre 5.5
is plotted to compare the relationship of Tm and size in the
simulation with experimental data in Eiggre 5.3, and gives

good agreement.

Energy E(red. unit)

C,(arb. unit)

2

I
.0
In
to

I
.°
0
as

-0.96

-0.98

- 500

100
50

10

95

 

ITIIIIIIIIIIII'IIrIIIII

I I I I I I I I I I I I I I I I I I I I I I I I I

N=5 381
i 177 72

A
D
v

lljllljllllllllllJJlJlLL

 

 

 

 

 

. 1
3 1

l l I I l l l I l l I l I I l I l l l I l l l I l
0 250 500 750 1000 1250

T(i()

Figure 4.4: Caloric curves for clusters N-SS, 177, 381 and 725.
ineachcurve,thetramitionpointisindicated(a)£v.s.1‘(b)c,,
v.s. T: The melting temperature is defined» the temperature at
which C. reaches maximum.

96

 

‘

[III

IIIIIllIllllllIllI'lIIllr

 

 

IIIJJJIIIIIIIIIlllllllllllll

 

-

 

10 20 30
Wk)

Figure 4.5: T.(r) from MD calculation: the four points are for the
melting temperatures 0fclusters N855, 177, 381, 725. The soihd
lineisatheoreticalfitfromequati0n(4.22).

A
O

97

An alternative criteria for the phase transition is to
look at the mean square displacement(m.s.d.) as a function
of simulation time t because atoms in the liquid show
adiffusive behavior while atoms in the solid do not. The

m.s.d. is defined as:

N
(w)2=%2[ri(t+7})—ri(7;)]2 (4~25)

ri(t) is position of atom i at time t. In calculating 5d,
the system reaches equilibrium first and the position of
each atom at time T1 is recorded as the initial condition
and results for an additional time t is used to calculate
(5d)2(t). t is fixed while Ti is arbitrarily chosen from
different values so that (Ed)2 is obtained from averaging
over various configurations. In our calculations, typical
data is averaged over 100 configurations and t is 10,000
time steps. The simulation is focused on the N=725 cluster
first. Eiggrg 4.6 gives (6d)2 as a function of t at various
temperatures. The picture shows that (5d)2 versus t has
different behavior at low temperatures from high
temperatures. At low temperature, (Ed)2 is small and
independent of simulation time t and increases monotonically
with temperatures; at sufficiently high temperatures, (8d)2
is of roughly linear dependence on t and 6d has magnitude of
order of distance between atoms. The small T behavior is due

to the oscillation of an individual atom near its

(60'

tad)”

98

 

0.10

A'I'IIVIVIIITIV'fTIVITIVI

 

0.4

 

 

 

 

1' VIII IIIT 11st ifiiilatirtltwtrlUVIV 111m 1111‘11

l I l l 1: l I 1
q

.E 1

(a): T-880K 034::- (b): T-880K '1

'll’ "

‘1- d

:11- .

u— -

an 4

db d

db 1

4h ..

u- -

qr- q

d- d

d1- «1

d)- 1

-‘* cad- -

7' fl. 1

d1- «I

d- «I

..IJJJIJAAJIJAHILJJJlJJAAIJA.

tiwrfirrrr‘rrvv‘rw 1‘1111*fi

(c): T- 1070K (d): T-lOBOK

llllLlLlllllllllllllllLll

 

been;

 

 

LI 1 l 14 A j l A 1 U1 1 1 AJ 1 .l l 1 ll 1 l l j l j. l l j l 1 1J .LJ 1 l 1 l 114 .l 1 ll
3 4 e a 10 0 s 4 o e 10
time step(x10’) time step(x10a

Figure 4.6: Mean square displacement (Iii)2 at various
temperatures for the cluster N=725. At low temperatures, (521)2
is small and independent of time. As temperature increases,
((21)2 increases in magnitiude and also has a linear dependence
on simulation time t, showing difl'usion behavior.

99

equilibrium position. The potential energy deviation 5V from
the ground state is approximately proportional to (5d)2 from
the harmonic approximation. The thermal average of kinetic
energy is roughly equal to that of 5V, therefore (Ed)2 is
approximately proportional to temperature T. The high
temperature properties of the cluster are a consequence of
atomic diffusion. High T allows structural changes which
create vacancies in the cluster and activates energetic
atoms to move into those vacancies. To quantify, the
diffusion constant is approximately measured from the slope
of (5d)2(t). The temperature dependence of D is shown in
Eigurg 5. . In the figure, D is zero or immeasurable at very
small T. There is a large change in D at the melting
temperature Tm. But it is also very remarkable that D has
non-zero values even at temperature well below Tm. This
implies the presence of diffusion in the solid-like phase, a
seeming contradiction with conventional wisdom. Tm for
cluster N=725 is roughly 0.0243(1070K) while the temperature
T. at which observable diffusion phenomena occurs is about
800K. This phenomena was not observed for clusters with
periodic boundaries and therefore, one has to look into the
detail of surface atoms. Surface atoms are less bonded and
easier to activate than bulk ones. We need to calculate mean
square displacements of bulk and surface atoms separately.

It is impossible to define a complete surface atom because

Math. unit)

100

 

 

 

 

-..1....l-s..l....l..fil....
.

b 1

, .
0.10— 11-725 -
. .

. .

. .

1- a
0.05— -
b d

. .

. .

D d
0.00-r 4r -
. .
....11-4.1....1....1....1.-..

 

500 000 700 1100 000 1000 1100
run

Figure 4.7: Difiusion constant D as a function of temperature T.
AtlowT, thedifi‘usionisabsentDw; thereexistsabig changein
D at T=Tn; it is also remarkable that D¢O for T4”.

101

it has a certain probability to be embedded in the bulk at
certain time. But it is unambiguous to define a group of
bulk atoms in the center of the cluster at low temperature
since the time needed for those atoms to move to the surface
is much longer than our simulation time. Therefore, in
defining the bulk atoms, the position of center of mass of a
cluster is calculated first and then the distance of each
atom to the center. A threshold value of distance is set to
identify the bulk and surface atoms. Atoms within the range
are called bulk atoms and those outside of the range are
surface atoms. Using the same technique, we repeated the
calculation of 5d but treated the bulk atoms and the surface
separately. The comparison is made at different
temperatures(£iggrg_gr§). At very low T( less than 700K), 6d
for the bulk atoms is similar to that of the rest and
independent of time t but has a lower value. At a medium
range of T (800K<T<1000K), (8d)2 for bulk is still
independent of t but (6d)2 for the rest exhibits diffusive
behavior, i.e. roughly proportional to t. It is clear now
that in this temperature range, diffusion does not exist for
bulk atoms but does appear to the rest of the atoms
including surface atoms. At high temperatures(T>1100K), the
bulk atoms and the rest do not act differently because then
it is meaningless to define the bulk atoms and surface atoms

as each atom diffuses around the cluster. Again from the

102

 

IIqq-I‘dd—dddd—dddd-dd‘dl

J

l

 

 

IIIIthIlIIIIlITTIIIIIII
I
(11)-r4001:
O

 

—
D
b

ITII‘jIII‘IITWiTIII‘jTTIl

IIIId-‘ddd_dqd‘—‘CII-‘dd d I‘dd—dd‘q—dzd-d‘qdd—qud-#dddl
1

.4

(d): 1400:

 

 

 

' l l l 1
IllJlellll 1111 1111 1111

all»
)-
I-
b
db
)-—

 

dr-
cil-
II-

fee“
all-
db
1-
an-
1
III-
alr-
eoaJ— r
h
db
d1-

TIME

 

 

ITIIITIIIlIIIIlIIII‘IjIII

 

thD—Pbtbbfibthbe_PPtl

 

1

0

10

(O'T-lm
lllllLllllllllllllll

4 I

timestepOtiO')

 

 

 

 

11

J

1

1

J

11

I

n 1

l

1. .1

J

m .

l

\...r 1
(w

J

4

1

j

.i

J

I

l

 

mum 4”“
a 2 2.63m. m m a

0.1

.95

i
1

 

 

'
0
._ ‘ I. ..
1m 1......
I 1
1 l
I l
1 1
i I.
.- 1 \l
r ._ I
o
I J 1
I m 1 x
.8 1 18W
. 1
I h L
4 s.» u
I
I ..t ...m
I l
I 1
T 1
T l
.2 L...
. h
I l
FEEL-LP.“
o

 

u w
1.3:.

Figure 4.8: Comparison of mean square displacements of surface

atoms(upper solid curves) and built atoms for the cluster N=725.
(ti!)2 ofthesurfaceatomshaslargervaluesthanthatofthe

bulk; The bulk atoms do not have difliisive behavior at T4 ,

andthesurfaceatomsaredifiusiveevenwhenthetemperatureis

well below melting.

DJnrb. unit)

103

 

0 a
. 10 F fi| ' I Y
N-725
10" , 0 °
0
10'2 ’
10'3 .

10"

‘0.5 1 A L A l L A A A l 1 A j_ 1 l A 1 1 J

 

 

0.010 0.010 0.020 0.025 0.000
‘I‘( orb. unit)
10° , it ,v 10.: V
l l l i I
l-‘Ifl
10" o
O
O
10" °
10" ,
10“

‘04 L 1+1 1 l 1 L l l A L l l 1 L 1 1 I 1 l l l I 14 L L
40 {I 50 I. O Q ‘N
t/‘K 0:5. unit)

Figure 4.9: Surface difi'usion constant D, as a function of T: D,
numekmhnfly ummnnes‘wfih I; Dy k nonwonkhhg vdum
1>800K; the lower plot is Ds v.s. l/T, showing the rough linear
dqxndwme.du:shperxmummsuubnmnkn untheawfiwnmn
emngy

104

relationship of (5d)2 and t, one is able to estimate the
diffusion constants D of the bulk atoms and D. of the
"surface" atoms. Eiggze_5&2 gives D, as a function of
temperature T. It is seen that the surface diffusion
constant D. decreases as temperature is decreased. Below
800K, surface diffusion is very small. The same calculations
were done on a cluster N=381 and similar phenomena is
observed. The simulation data are summarized in Eiggzg_gélg.
It is also found that the temperature below which no surface
diffusion is observed is roughly equal to 800K, the same
value estimated from the pervious cluster calculation. But
when applying the calculation to N855 cluster, the surface
diffusion and the bulk diffusion seem to take place
synchronously and no diffusion occurs prior to melting. For
the intermediate size cluster N-l77, we did observe surface
diffusion at the temperature below but very close to the
melting temperature. Because of the difficulty in separating
surface and bulk atoms at temperatures near Tm, our
calculation results are not conclusive. In summary, surface
diffusion exists in large clusters( N>200) at temperatures
below melting. But if the temperature is too low(T<800K), it
is very hard to observe. T.(about 800K) does not change much
for various size clusters but the melting temperature Tm
decreases rapidly as N decreases. As a result, there exists

a critical cluster size NC and when N<Nc the melting

O.‘

8

(«P

0.1

(60'

0.4

1.6

(MI'

0.0

105

 

 

 

 

 

 

 

 

 

 

 

 

111W Yj‘Ir 111W I‘IWT 1T1! Tifr the test its III’T 11
p P
_ I I I I i I I T‘ I I
I- “l: '1
y. Ii «I
_ 0 ° —
. (I). 'r-eeox 0:01? (I). r-eeox ‘
# db 1
p . ‘
r. ..L' .
— 09- -
b db 1
I- q
E a- -
h «I
— 0am- -
p h 1
p j? I!
P ‘P 1
p Z ‘P ‘
- 0m— -
p ‘D ‘
P 1. 1
t t I ‘
. jj j l LLJ l l l 1 All 114144 lj lLJfi
' Ijj I I I I I‘I I IW‘I I rI‘I I I I I I

 

 

 

:

(arr-070: (an-0201: {

/%

J -:

4,

(a): r-oeox (on-10101: j

;

1

1 —I

AIAALLJIAIAALLILAAAIIJLIl lJlllJlllllllllllJllllll A:
O 8 O I 10 0 8 4 I I 10

‘
tune ”(x105

time mp(xro')

Figure 4.10: Mean square displacement calculation for the cluster
=38l: In the calculation, the bulk and surface atoms are treated
separately. This cluster has similar difl'usion properties to the

cluster =725.

106

temperature of the cluster is too small (T<800K) to have
surface diffusion below melting. A very rough estimate of

the critical cluster size gives Ncslso.

WWW

Current researches on the dynamics of metal clusters are
limited to the calculation of vibrational frequencies by
utilizing fast Fourier Transform(FFT). The direct structural
changes a small metal cluster may undergo have not been
studied. On the other hand, the advanced microscopes enable
experimentalists to observe single atomic motion(Ajayan and
Marks, 1988,1989). In this section, an effort is made here
to narrow the gap between theory and experiment through
investigating the shape change kinetics from a non-
equilibrium state to an equilibrium state.

Simply, a system in non-equilibrium means that the
system deviates from an equilibrium state, which possesses
the lowest possible energy. In reality, a system could have
many states with energies very close to the ground state.
Those states usually have extremely large energy barriers
between them. At finite temperature, the free energy F can
be roughly considered as a function of configurations. F has
numerous minima and each of them represents a meta-stable
state. Large energy barriers separate those states and makes

it difficult to transfer from one to another. The higher the

107

temperature is, the more likely it is for the transition
from a metastable state to another to occur.

We study small clusters( sto ) first. The initial
configurations of clusters are set to be of an elongated
shape. The goal of our calculations is to observe the shape
change of the cluster. Initially, any portion of a cluster
will have same temperature. It is realizable when the
temperature is small because the typical number of time
steps necessary for a system to complete heat exchange is
much smaller than the time entailed to alter its shape
significantly( as will be seen below). We concentrate on
cluster N=41, which initially is approximately a 4 by 1
cubic rod(£igu;g 5,;1). We adopted the steep quenching
method(Stillinger and Weber, 1982) to freeze the positions
of the active corner and edge atoms in order to maintain the
initial elongated shape. In the quenching process, we reset
the velocities of all atoms to be zero after a short period
of simulation. Consequently, the cluster quickly reaches a
very low energy state and keeps its elongated shape. The
shape parameter is defined by the radius of gyration r

9!
written as:

l
r82=XI-Z[r,(r)-rc]2 (4.26)

where rc is the center of mass (c.o.m) of the system that is
stationary because of zero total momentum; ri(t) is the

position of atom i. Geometrically, rg measures the position

1.08

 

 

N-41.r'.-1.43

.'

 

Figure 4.”: Typical initial and final configurations of cluster
N=4| in studying the kinetics of non-equilibrium shape change.
The upper figure shows an elongated structure which has a value
of r =1.65; the lower figure is the configuration afier reaching
equiiibrium with rg=l.43.

109

distribution of atoms with respect to c.o.m, therefore an
elongated cluster has a larger value of rg than a round and
compact one. With the restriction of constant volume, it
also measures effective surface area. A configuration of
large rg is unstable and tends to decay to states of smaller
rg values. Our MD simulation is to look at the time
dependence of r9. £1gg;g_1;12 presents the time change of rg
for a cluster N=41 at temperatures T=300K and T=360K. To
interpret the phenomena explicitly, two time scales tc and
tq are defined as indicated in the figure. tc called the
completion time,is the time step at which rg is
approximately at middle of the initial and final value. tq,
the decaying time characterizes the time step needed to
finish the transition of rg from high to low. In Eigure
5‘12, tq is of order of tens of thousands of time steps
while tc is roughly one million time steps. tc is strongly
dependent on the temperature. At very low temperature, an
elongated configuration in a metastable state may not alter
its shape for such a long time that a time scale loans has
no impact on the overall profile of the cluster. I did
simulations on an N=41 cluster at about half of its melting
temperature(24OK), and the simulation was extended to 107

time steps. No big change in r occured. It is found that tc

9

also has a dependence on its initial configuration. We fix

temperature and change the initial condition for the cluster

110

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

107 I I I I l I I I I l I I I I I I I I I I I I I I q
N-41.'r/'r.-0.05 j
R“ _"
t 1
1.6 — . —-i
I- d
b
t.‘ I. l l l l l l 1 J l l l l l l l 1 l l l 1 7 1 1
0 000 400 000 000 1000
10? .. I I I I I I I W I I I I I I I I I I T .1
: N-41.t/1‘.-0.00 j
r I
1.0 - —
I .
1.0 — -
P u
1 ‘ P l l l 1 J 1 l 1 1 l l L J j l l 1 l l
' 0 100 200 300 400

time step(x10')

Figure4.:12 Radius ofgyr'ationrg(t)asafunction oftimeat
T=300KandT=340K for N=41. Ing,thesetwofigures r g(t)has no
significantchangeuntilalargetimescalet asdefinedinthe
upper figure, but completes compactification within a very short
period tq.

Pu.)

0.25

111

 

0.20

0.15

0.10

0.05

 

0.00

 

\

JillllLlllllJllllllllllL

l l llllJll I l llllll

 

 

 

10

20 50 100 200 500 1000 2000
11.0: 1000)

Figure 4.13: The distribution of tc for N832 clusters. We change
the initial configurations and fixed the temperature. The statistics
is done on 100 configurations. tc has a very broad distribution;
the figure shows a linear-log plot which yields a reasonable
statistics.

112

N-32 and find that tc has a distribution so broad that a
simple statistical analysis failed. Eign;g_1L11 shows a log-
linear plot of distribution of tc, indicating that only
logtc has reasonable statistics. As to the system size
effect, we fixed temperature and chose clusters N=59 and 67
for comparison with the cluster N=41. We did observe that tc
increases as N increases. Nevertheless, the melting
temperature is also increasing with N, so the temperature
difference from melting point is also increasing with N.

In large clusters, it is more difficult to observe
substantial alteration in rg at low temperature, therefore
we restrict our calculations to temperatures close to
melting. The clusters chosen are N=711, 613, 515, 417, 319
and 213. First we estimate the melting temperatures of these
clusters and find that Tm=0.0242, 0.0238, 0.023, 0.022,
0.021, 0.019 respectively. Again, the starting profiles of
clusters are set to be elongated and the radius of gyration
is calculated as a function of temperature. We run the
simulation up to 3,000,000 time steps. Our calculation
results are presented in £1gu;g_1;15, which shows how rg
behaves at different temperatures and system sizes. For
cluster N=613, the temperatures are chosen to be T/Tm=0.96,
0.92, 0.86. At T=0.96Tm, rg has the abrupt decay as
described in the previous paragraph. But at T=0.92Tm and

0.86Tm, rg decays rather smoothly and gradually. It is also

113

noticed that the decay at 0.86Th is significantly slower
than at T=0.90Tm. For the cluster N-SIS, the temperature

chosen are T/Tm=0.90,0.85 and 0.78. At T=0.90Tm, r has

9
arrupt decaying behavior while at T80.85Th and 0.78Tm, the
decay is gradual. For the cluster N=711, at temperatures
T/Tm=O.92,0.88,0.82, the decay of rg is always gradual, but
we expect as the temperature is close enough to the melting
temperature the abrupt decay behavior will eventually show
up. For the cluster N=213, we did not observe the slow decay
behavior. If temperature is close to melting, the abrupt
decay is observed and if temperature is low, the overall
shape is hardly changed. Accordingly, two distinct phases
can be categorized. Let Phase I label the gradual decay of
r9, and Phase II stand for the abrupt decay. Then we
observed that at temperatures very close to melting, the
cluster is in Phase I and crosses over to Phase II as
temperatures decreases. Because the slow and gradual decay
behavior of rg exists only in large size clusters, it

implies that the decay mechanism is surface diffusion below

Tm, which is absent in the small clusters.

8.1

B i 2 8 8 22!:

’3

114

 

l J I l l l l l J I l l l L I 1 l j I
I I I I

I I 1 I j fi I. I I T I I I I I I I I I I I I I I
(am-210; r/‘r.-0.0e

um Mario)

1 l l l l l I l J l I
I

(coat-eta; r/r.-0.00

l
ar'fluo'H'uo 1..

(om-210: 'r/‘r.-0.00

1 1 J l 1 l 1 1 I 1 1 J 1 l 1 1 1 1 l 4 1 1 1 ,
400 000 000 1000
time Itep(x 100)

115

 

4.: (b):N-515: T/‘l'.-0.?l

4.0

l I I I I J I I J I
I

u l l 1 1 1 I
I I I I I I I I I I I I I I‘I

4.0 (um-51o: rn.-o.oo

 

 

 

 

u mac-515: 'r/r.-o.oo

8.0

l 1 1 L .l 1 l l l l I 1 1
ooo ooo zooo

&‘ 7 I I I I l I I I I l
coo .
tint utop(x600)

116

 

 

“ I I I I I I I I I I I I I I If I I 1 I V ' j I
u (om-us; 'r/‘r_-o.oo
4.8
‘1
I : : : l I I I I I I I I I l I I I I I I I I I
‘ * I I I I w I I I I * I I I I

(ow-61:; 'r/‘r.-o.oz

 

 

I

d

d

.

db d
: llllllll IIIIIIIIIIIIIILL‘E

 

an (ow-us; 'r/‘r.-o.oo

ooo ”0 10M

117

 

8$££§8$ t t C_’." I: t 3‘: '3

II I I I I I I I I I I I I I I I I I I I I I I I I I I I I

(aw-711; r/r_-o.oa

:...

mar-711: 177,-one

1 . ' 1- ..
(0:11-71 1; 'r/‘r.-o.oo

I I I I I I I L l I I I I J I I I I I I I J I l
ooo won I.” aooo aooo
tune otep(x1000)

Figure 4.14: r‘(t) at difierent temperatures T<Tm for the large
clusters N=213, 515, 613 and 711: the temperatures are
normalized to their melting temperatures. r‘(t) decays distinctly at
'1' close to melting and at T a little away from it. When T is very
close to Tm, r‘(t) decays abruptly compared with a slow and
gradual decay at lower temperature.

L I

I

 

118

So far, only clusters of solid—like phase have been
considered. When the temperature is raised to surpass the
melting temperature, atoms in the cluster start to diffuse
and therefore, the transition from a non-equilibrium state
to equilibrium becomes much quicker than that in a solid
phase. Since the time scale of the transition is compatible
to that of thermal equilibration, the term isothermal shape
change is inadequate here to describe the prompt
equilibrating process. However, it is still feasible to
utilize rg as the parameter to look into the shape change
because the shape alteration of clusters will certainly not
change the phase of the cluster, in other words, the cluster
will remain in the liquid phase. A single configuration of
rg as a function of simulation time for a liquid cluster has
large fluctuations as seen in Eigg;g_5‘1§a. But it is clear
that the statistical average yields a rather smooth curve
which has a well defined behavior(£igg;g_5;1§b), showing
roughly exponentially decay(£1gg;g_1&;§g). From the averaged
decay curve, we measured the decay constant I. We measured 1
as a function of temperature for three different sizes of
clusters(N=4l, N=213, and N=613). The results are plotted in
Eigg;§_5&1§, from which two conclusions can be drawn: (1) 1
decreases as temperature increases, and it reaches a stable

value at high temperature. (2) The stabilized value of t is

119

 

1.7

 

“-
1.5
1.1
m
8.7.
“-

 

IIIIIIII'IIIIIIIIIIIIII‘ITU'fi

(gm-11; r/r.-1.e

III IIIIIII
109111101011

8‘
p-
—b
d

mar-21:: r/r.-1.o

IIIIILIIIIIIIIII

 

 

L”

0.10

120

 

II"'

'I

I

 

 

IIJLI

Oahu-41: 177,-1.5

 

 

Figure 4.15: r'(t) at 151'“, for the clusters 14-41213. r t)
deaysapproxirmtelyapomtiallybutwithbigfluematiou.
avengebehaviorofrgtuasawelldefineddeenytimeoonstamt
whiehisdefinedutheslopeofthelog-linareurveabove.

121

 

 

 

 

 

1.1. III. 11.! .111 1‘1. 1“
_ _ — _ _ . 1
I I
I
I I
I I
I II
I I
I
I I
I I
II I“
I . II
_. I
t j'
I I
I In“
I I
'
I 1 1
I 4 II
. fl .
d. .
I. ( II
I I
I.
I I
I II
1 I
I It”
I I
U
I I
I II
. I
beb—tpD—tbb-IIFD-bbe-Dbbb
. . I . ' 1.

A972 5b

a?

d-q‘.‘-IICI-CqCI-IIII

(IN-818

5373.57

_

 

(GIG-QI‘I-q‘Iq-dd‘I-CI

 

I
I
I
I
II“
IO
II
III
I
I
I
I
I
In
I
I
l u
.. a
.. M
O
(

 

 

I“

l
1000

 

decreaseswith‘l'fitmblinsatsuficientlyhigh'r. Thestablized

Figure 4.16: 1 mm for the clusters N-41, 213, 613. 1
value is a weakly increasing function of system size.

 

122

a weakly increasing function of system size N. The
conclusion (1) is consistent with the phenomena that the
equilibration constant of a liquid droplet is weakly
dependent on temperature.

In summary, the molecular dynamics simulation is carried
out to study the kinetics of non-equilibrium shape change in
clusters. We utilized the radius of gyration rb as a shape
parameter to monitor the equilibration process. It is found
that rg behaves distinctly in three temperature regions for
large clusters. Above melting, rg decays roughly
exponentially. Near and below melting point, clusters
undergoes rapid compactification at some characteristic
time. In large clusters at low temperatures, rg decays
rather smoothly and algebraically slowly. In small clusters,
rg rapidly compactified at temperature below melting. Those

decay behaviors of r are summerized in Eigggg 5,17 that

9

shows how rg decays at various temperatures and cluster

sizes.

'1‘(K)

1200 - . . I . . I j j I . 1

r

_ Phase I
1000 - ..

h d
800 [- Phase II -

. Phase III .
600 - '1
I l I I I I I l 1* I I

400 ‘ ‘

0 200 400 000

N

123

 

 

 

Figure 4.17: Dynamic behavior ofgold clusters. In the figure,
threedistinctbehavionofdecayofrgthreidentifiedlnl’lmel.
itoorrerpondstothecueofliquidcoaleeonee; Pharenmnds
fortheabruptdecnyhehafiorlndl’hasemisslowdecaydueto

surhoe difi'usion.

 

CHAPTER FIVE
QUANTUM TRANSPORT ON 2D PERCOLATING MEDIA

Wises.
5.1.1 Anderson localisation

The electronic states in a perfect crystal material are
extended with Block wavevector k. An incidental defect in
the crystal can be viewed as a scattering problem of
electrons on the defect. When there are many defects
randomly located in the crystal, an extended electronic
state may be altered profoundly in nature( Anderson, P.W.,
1958) by the disorder and collective effect of defects.

Roughly speaking, there exist two broad distinct classes
of disordered systems: i)systems with substantial, cellular
or compositional disorder and ii)systems with positional,
topological or structural disorder. The first class of
disordered systems include the famous Anderson model
(Anderson,P.W., 1958,1961) which uses the following

Hamiltonian:
H = 28,0;0," + ZJM.(a;a,, +a;am) (5 . 1)
or sun!

where a, is the energy level of site m, and J”, is the

overlap or hopping integral between site m and its nearest

125

neighbor m'. The disorder was introduced by making the site
energy be uniformly distributed over a certain range, -W/2
to W/2. J“, was assumed constant. Let V stand for the energy
bandwidth in the absence of disorder, and the ratio W/v
represents the magnitude of disorder. When W/V is
sufficiently large, all electron states are localized,
having the asymptotic form(Lee etc.,1985):

|W(r)l~eXP(-lr-ro|/€) (5-2)
where § is the localization length. When W/V is small, the
localized states show up only at the bottoms and tops of an
energy band. In an energy band, there exist energies called
mobility edges separating localized states and extended
states. As energy approaches a mobility edge, é diverges:

545-1%” (5-3)
where the index v reveals the critical behavior of E and Ec
is one of those mobility edges.

In terms of Hamiltonian (5.1), an alternative way of
introducing the randomness is to let a, be constant, and J”,
random, which corresponds to the situation of off-diagonal
disorder and the second class of disordered system. The
simplest model with off-diagonal disorder is that of quantum
percolation which takes p as the probability that a bond
exists or 1-p as a bond is absent. The raised question is
whether the localization caused from off-diagonal disorder

falls in the same universality class as the well studied

126

site disordered case. To gain a quantitative understanding
of the question, one obvious approach is to carry out
detailed numerical calculations or simulations and compare
the results with the theory. In this chapter, we performed a
detailed calculation and analysis of conductance of two
dimensional percolation systems. Two key issues are
addressed in the chapter: (1)all electronic states are
localized in 2-d percolation systems; (2)the wavefunctions
of those states might be superlocalized near classical
threshold. In the rest of the section, I will go over some
fundamental concepts in percolation theory and localization
theory. In §§gtigg_§&z, a description of two methods of
calculating conductance will be outlined. Then in Section
5&1, the numerical calculation results are presented.
Finally, conclusions and open questions are discussed.
5.1.2 Scaling of localization:

The nature of Anderson localization was unclear until
the scaling theory of localization was discovered(Thouless,
1977). The divergence of the localization length 5 near Ec
implies one could employ the renormalization group used in
analyzing thermal phase transitions where the correlation
length diverges as the temperature approaches the critical
temperature. When applying the real-space renormalization-

group method to localization problems, one divides the

127

infinite lattice into finite cells of bd and thus the
Hamiltonian is written as:

H=H0+H, (5-4)
where Ho includes all site-diagonal terms and off-diagonal
terms which connect sites m, m' within the same cell. The
remaining off-diagonal terms are included in H1. The
eigenvalues E; and normalized eigenfunctions PD“) of Ho can
be found by exact diagonalization. The objective of
renormalization is to find a mapping Hamiltonian

characterized by E;, the renormalized local energy and Jun

the renormalized nearest neighbor hopping cOnstant, i.e. the
relationship between (EMJM) and (£_,JM). From each cell,
one picks out one state which is included in the subspace D
of Hilbert space that Ho spans. Introducing the projection
operator PD of I‘I’) onto D, I‘PD)=PD|‘P) Where PD can be

expressed as:
PD =Z|¢a)(¢a| (5.5)
050

Defining the Green's Function

G=(E-Ho)"(1-PD) (5.6)
then

|‘P)=|\PD)+GH,|‘P) (5.7)
Defining another operator 0 by l‘P) = QI‘PD), and then VD by
VD = H,Q, one finds Q=1+GHIQ and therefore

VD=H,+H,GVD (5-8)

128

Noticing VDI‘I’D)=H.|‘I’), it is found that
(H.+VD)|*PD>=EIWD> (5.9)
One can obtain Vb by expansion in series of H1. Keeping only

the linear term i.e. VD-Hl and defining HD=H°+VD, we have

H=ZEaa;ad+ZJw(a;ad+a;aa) (5.10)

and

where (a,a') are nearest neighbor cell indices and Jaa.=(alHl

la'). Principally, one could find the recursive

relationships between (Jaa.,Ea) and (Jm.,em). Consider the

diagonal disorder situation, one could make the

approximation that J;,=.Q, are constant and.Eh has a new

distribution whose width or effective disorder is written as
"2
W‘:(<E§)-(E)2) (1., (5.11)

and the energy E is changed to E':

E=(E-(Ea))/J,J (5.12)
Then the scaling for localization length reads

€= 6(W'.E')=b"§(W.E) (5.13)
If one considers the case that the energy lies at the
center of the band, according to the symmetry, E=E'=0. In
other words, in the (W,E) plane the E=0 axis is invariant
under RG transformation. Thus the recursive equations become
W'=f(W) alone. The fixed points are given by solving the
equation W=f(W). Obviously W=O is a stable fixed point
because the renormalization process should not introduce

disorder in to a perfect system. The non-trivial solution

129

w=w¢:o, if it exists, represents a critical degree of
randomness above which all electron states are localized. In
that case, W - We is an unstable fixed point at which the
localization length a diverges.

More quantitative analysis of scaling of localization is
to consider the behavior of conductance F or reduced
conductance g defined as g=l" /(e2/h) as a function of system
size L. When randomness is small, the electronic wave
function is extended and is nearly plane wave-like. The
conventional transport theory is valid and the ohm's law
gives:

g(L)=ai."‘2/(e2/h) (5.14)
here a is the conductivity viewed as a local quantity
independent of system size L. Equation (5.14) is the ohmic
limit of g(L). In the limit of strong localization, the
conductance is due to hopping from one occupied state to an
unoccupied state. The conductance is then proportional to
the probability that an electron hops between two
energetically close states which have a typical separation
L. This probability is proportional to the overlap of two
localized wavefunctions, thus in the strong disorder limit,
g(L) has the scaling form:

g(L)=gOCXP(-L/§) (5.15)
go is some microscopic conductance value that is of order

one .

130

For a particular degree of disorder, g(L) evolves
smoothly as L increases, crossing over to either an ohmic
region or to a strong localization limit. Here we use the
argument given by the gang of four(Abrahams,E., etc. 1979)
that assumed the quantity

fl=dlogg/dlogL (5.16)
is a function of g alone, not explicitly dependent on system
size L. When 9 is very large, equation (5.14) is used, then
B is easily computed: B(g)=d-2. Particularly for two
dimensional systems B=0. If g<<go, using the strong
localization limit(Equation (5.15)), B-ln(g/go), which still
depends on g only. Perturbation theory in the large g limit
proves that fi=d-2-a/g(Bergmann,1984) where a is a
positive constant. As it is presumed that B is a
monotonically increasing function of g, B is always negative
in one and two dimensional disordered systems. The non-
positivity of B implies that any finite degree of disorder
will cause non-conductivity in 1d or 2d macro-systems.
However, the conclusion is still controversial(Azbel,1991;
Soukoulis and Grest,1991; Stein and Krey,1987) in the
situation of two dimensional disordered system where the
localization length is very large at small disorder(Lee
etc.,1985).

5.1.3 Scaling behavior of percolating systems:

131

The most important feature of the percolation model is
the occurrence of an infinite cluster(Stauffer,1985) of
sites(or bonds) for p>pc, where pc denotes the percolation
threshold value. P(p) is an order parameter defined as the
probability that an arbitrary site(bond) belongs to the
infinite cluster. When p<pc no infinite cluster exists, thus
P(p)=0. As p is reduced to approach pc from above, P(p) has
the scaling behavior:

P(p)~(p-p.)‘ (5.17)
indicated by exponent B. The correlation length and
classical resistivity both diverge at pc.

§,~(p-P.)"’ (5.18a)

p.=l/a.=(p-p.)“‘ (5.1810)
Near pc, the structure of the infinite cluster is fractal
with dimension d . d, can be obtained from a simple scaling
argument. Consider a block of size L in the infinite cluster
and L satisfies: a<<L<<§p, where a is lattice constant. The
weight(mass) of the block is proportional to P(p)L" = Ld’ . Let
L==§,, in which case we have the following relationship:

(p-p.)”~§i"" ~(p-p.)"“""” (5.19)
and therefore d, =d-fl/ v.

Diffusion on the infinite incipient cluster at pc is
anomalous and the mean square distance of a random walker

from the origin after time t behaves as 03>~t”“ (when dw=2,

we get normal diffusion.) The diffusion constant on the

132

cluster is defined as D=§d<r2>/dt and D relates to the dc

conductivity(Harris,1984) by
a=eanlkT (5.20)
where n, the density of charged carriers, is proportional to

P(p). Thus, the simple scaling argument yields
D~(p-pc)""~§: with 6=(p—fl)/ V. In general, for a scale L<E,,

D scales as D~§:f(L/§,). When L is small, D should not
depend on :9, which implies the function fo)~x’ as x—9O.
Therefore D~dr2/dt~r’9 or

r’(t)~12"”’ (5.21)
This expression is useful in discussing of super-
localization( Levy and Souillard,1987; Brooks and Aharony,
1987) of electron states on a fractal lattice. The topic

will be briefly discussed in last section of this chapter.

MAW

5.2.1 Linear response theory:

Many experiments in condensed matter physics measure the
linear response to small external fields such as mechanical
forces, magnetic fields, electric fields or optic fields on
materials. In electrical conduction, a time dependent
external electrical field

E=Eoexp(iq-r-iwl) (5.22)

133

is applied and the induced current is measured. Linear
response suggests that the induced current is proportional
to the applied field:

J(q.w)=0(q.w)-E(q.w) (5.23)
where o is a second rank tensor and the expression can be
viewed as the fundamental definition of the microscopic
conductivity. The dc conductivity is obtained from taking
the limits qe+o and m~+0 in that order.

When a field is applied, the Hamiltonian consists of two
parts: one is the unperturbed term Ho and the interaction
term H' is expressed as:

H'=-%Jd3rj(r)-A(r,t) (5.24)

lA(r,z)=-iE(r,z) (5.24a)
C a)

The current density relates to average velocities through

J(r,t)=eZ(v,) and v, =-};[p,—§A(r,)]. Therefore, the current

density is sum of two parts: J231+Jz, with J} and 55 having

the following expressions:

 

. e
1. =(:(r,t))=;Z(p.-) (5.25)
J1 is the displacement current and proportional to B
- 2
J, =m°e E(r,t) (5.26a)
:nw

The expectation value of j(r,t) is expressed in the

interaction representation

J2(I'.t)=('//'|j(l'.t)lV") (5.26b)

134

Iwfi represents wavefunction of the system at time t and it

evolves with time as:

IV)=S(-°°J)i W) (5.27)
lué denotes the wavefunction of the system when no field is
applied. For very small external fields, the 8 matrix is

expanded to the linear term of H', yielding

S(—oo,:)=[l-i/njdrH'(r)]+--. (5.2s)

Then the current J; can be expressed as a commutator:

Jz(r,t)=-i/hjdt'(w|[j(r,t),I-I'(t')|yl) (5.29)

Expressing H' in terms of the electrical field and the

vector potential, one obtains

Jz(r, t) = J; Idt'e“"”( M[j(r,t),j(q, t)|w>E(r,t)exp(-iq ~r) (5 . 3 0)

Now one has the expression for each conductivity component

in terms of the current-current commutator

0a,,(q, w) = % jdt'e'w-r)< ”I“; (q, t), j ”(q, t' )]| W) 1 no: 1'60, (5 . 3 1)

m

 

This equation is the Kubo formula for conductivity which is
fundamental in calculating the conductivity or conductance.
It is desirable to have the expression for conductivity
or conductance in terms of eigenvalues and eigenfunctions of
the Hamiltonian which describes the disordered system(
Fisher,D.S. and Lee, 1981). If the system is isotropic,

¢z¢==afiw and the real part of conductivity can be written

135

Re<a>=i-jdte'“(wu<q.t),j(q.0)1 w) (5.32)

The expectation value of the current-current commutator is
expressed in terms of eigenfunctions of H0 at finite
temperature T (T=0 is obtained by taking the limit

fl=llkT —>oo) . The expectation of the commutator is

lz_(e - -<)

where MD and h» are eigenfunctions of HQ with eigenenergies

)|e “- '~“" (5.33)

        

Eh and En respectively. Z is partition function given by
Z=ng (5.34)

Integration over t in equation (5.32) gives

Re(a)=_a:fz.(1—e-”°)Ze”- (m(j(q)|n)|’a(w-E, +E,,) (5.35)

 

Now we consider the disorder system of size L with two ends

attached to two perfect leads. The current flows along the

z-direction
J(z)=Id"'xj(x) (5.36)
j(p) is the Fourier component of j(x) thus as q—M)
l L
j(O)=L-7—,2Idz.l(z) (5.37)

Finally one takes the limit of B-MD

n>2

            

 

Re(a)=-a—:—dIdE'z 0‘ 6(E'-E,)5(E'+w—E,,) (5.33)

The dc conductivity is attained by taking the limit of w=0.

This expression directly relates the conductance to the

136

eigenvalues and eigenfunctions and is our starting point in
doing numerical calculations in segtign_§‘;.
5.2.2 Landauer's formula and transfer matrices:

The calculation of conductance in one dimension is
greatly simplified by the formula which expresses the
conductance of a system of random elastic scatterers in
terms of scattering properties at a fixed energy. The
formula is called Landauer's formula(Anderson,P.W. etc.,

1980; Azbe1,1980) written as:
e2 T
=33; (5.39)
where R and T are reflection and transmission coefficients
respectively. Therefore if one knows the scattering matrix,
the conductance can be easily computed. The formula can be
derived rigorously in one dimension, however, the extension
to higher dimensions seems unclear. For multiple-channel
scattering, the generalized Landauer formula proposed by

Anderson etc.(Anderson,P.W.etc.,1980; Buttiker,1986) takes

the form

2
2 2|th
e a

 

r: 5.40
ZMI-ZIIUF ( )
1
When the channel number n is large, tif~1/n2, thus
2
e
r wz—Trtt’ 5.41
mo Znh ( ) ( )

This expression is accurate for strong scattering.

Therefore, it may be advantageous to compute the transfer

137

matrix first and then conductance can be fairly easily
obtained. As shown below, the transfer matrix can be
calculated column by column so this algorithm scales as L
while the matrix size of the Green's function method is L2.
Now let's consider in more detail about the propagation

of the wave function. We still consider the disordered
region to be attached with two perfect leads of infinite
length. Let ya be the wavefunction at the i-th column of the
disordered region. The tight binding model has the following
recursion relation:

Hy. =V.w.-.. +V.-.w.-. (5.42)
where H; is the diagonal matrix of column 1 and Vi is the
connection matrix between column 1 and i-l. This is

equivalent to:

-l _ 4
[w..]=[V. H. V. V51 v.) (5.43)
V) I 0 WM

The above expression gives the relationship between two
consecutive column wavefunctions. Then the transfer matrix T

is simply expressed as:

-l _ 4
“(V-1H. V.OV.-.] (5,44,

In order to get conductance, the transfer matrix has to be
converted into a scattering matrix S(Anderson,P.W.,1981).
The conversion utilizes the relationship: S=U"TU where U
is a unitary matrix and is determined by the Fermi energy

level.

 

138

W

In the section, I will use the Kubo's formula(Equation
5.38) for conductance to calculate the conductance of a 2-d
diluted square lattice which is described by the

Hamiltonian:

— x + y * +
H - 2V2“).u.ka)i + Vfia’Mafl + c.c.+ZEfiaflafi (4 . 45)
1* 1.:

where j and k are position indeces in the x and y directions
respectively and Ejk are site energies. In the bond
percolation model, V;,V; are hopping constants which are one
when a bond exists and 0 otherwise and Ejk1are constant. In
site percolation case, Ejk(are either infinite for absent
sites and zero for present ones. For a finite system of

scale L, the Hamiltonian is a szl? matrix expressed as

following:
(11“ H” H”)
H: If“ H.” I" H.“ (5.45)
\Hu 1L1!”2 Ha}

 

 

where each Hmn is an LxL matrix with Hhm,.Hmflmfl being non-

zero only if the nearest neighbor hopping is considered. We
start with the expression for the conductance of an LxL

lattice at frequency w. Equation (5.38) becomes:

2

n I O
G(m)=;Z-2-JldEz 6(E+m-EB)8(E‘-Ea) (5.47)

a..B

Z<alJ<j>IB>

J

 

 

 

139

where EWEB are the eigenvalues of H with 15°I <E<Ea where E

is the Fermi energy and |a),|B) are the corresponding
eigenvectors. Following P.A.Lee etc. (Lee and Fisher,D.,
1982) , one introduces the Green's functions 0* =(E-Hii17)"

where n is a small positive number. The imaginary part of

the Green's Function is G=(G’—G*)/2i. E and la) are

eigenvalues and eigenvectors of H respectively:

 

Hla)=Ea|a) (5.48)
and in terms of the Green's functions:

|a)=G*(E-Eaiin)|a) (5.49)
Taking the imaginary parts of both sides gives:

G|a>=(E-E:I)2+1fla> (5.50)

As co—)0, the continuity of current requires J (j) ought to be
independent of j, thus summation over j is trivially done,
yielding L2|(a|J(j)|fl)|2. Noticing that

. 17 _ _
m(£-E,)’+If-flaE E) (5.51)

 

One employs the stepfunction f(x) to eliminate the
restriction: E, < E'< E,

11(90):;16I‘1E0 f(E-E.+(D)f(E"E)

co g(alJGIBXBlJGIa)
=Z<a|JGJGia) (5.52)

 

The current operator J(j) of jth column is expressed as:

J(j)=ieZVl,(a;,a}_u—c.c) (5.53)
k

140

 

Em: A 32x32 percolation square lattice embedded in a
perfectsquarelattieebeltofwith32andinfinitelengthwith
p=0.6 in the percolation structure. All the data in Chapter five
are based on the model shown in the picture. The adjustable
pamneterssrepsndsystemsizeL.

141

one views the current operator to be a matrix of 2Lx2L(Lee

and Fisher,l982):

J(j)=[ .0 . ier] (5.54)
4er 0
where V7 is the vortex diagonal matrix with.(Vf»,=P:6;. The
conductance at m=0 is simplified to:

r(0)=e2/nTr[G( j, j)V;G(j-1, j-1)V;+G( j—l, j-1)V;G( j, 1')ij

-G(J'.j -1)V,‘G(j.j -1)V,‘ - GU -1.J')V,-‘G(j - 1.1)V,’] (5 - 55)
where the trace is over the site index k. The boundary
condition is set to be periodic along the y-axis and the two
ends x=0,x=L are attached to two semi-infinite perfect leads
(Eiggzg_§;1). From equation (5.55), to calculate the
conductance, one only needs to know G(j,j), G(j-1,j-1) and
G(j,j-1). To avoid a huge matrix representation the
recursive method is utilized. We start with a perfect
lattice with width L and semi-infinite length. We define the
left Green's Function (V(j)2(fl,(flj)'with the columns of j'>j
absent. G'( j) has the following matrix recursion
relationship:

G'(j)" = (30(1)l -Vj"G’(j—1)Vf (5.55)
where G°(j) denotes the Green's function for the isolated
jth column. Similarly, one could define the right Green's
function Gr(j) with all the columns of j'<j absent. The

Green's function of the jth column which connects the (j-l)

142

th column from left and the (j+1)th column from the right is
written as:
0031')" = (30(1)‘I - V,i.G'(j +1)V,.. - Vf..G'(j -1)V,‘.. (5 - 57)
and G(j-1,j) can be obtained from G(j-l,j-1) through:
G(j-l.j)=G(j-l.j-1)V,*G'(j) (5.58)
First, 61(0), the left Green's function for a semi-
infinite perfect lattice is computed. To find 61(0), one
just needs to solve the matrix equation
G’(0)‘I = G°(0)‘l -G’(0) (5. 59)
because Vf turns out to be a unit matrix. Because 60(0) and
61(0) have the same matrix form, the transfer matrix T which
leads to the diagnalization of 60(0) happens to diagonalize
the matrix 61(0) also. The eigenvalues and corresponding
eigenvectors for 60(0) are easily obtained and so that one
has the eigenvalues of 61(0) from the simple second order

equation.

1 l

where ’1‘: denotes the eigenvalues of G°(0) . 61(0)is
calculated from G'(O)=T"AT, where A is the eigenvalue
matrix of 61(0).

Once 61(0) has been calculated, one can build the
disordered region onto the semi-infinite perfect lattice
according to equation (5.56). After the L-th column has been

attached, another semi-infinite perfect lattice is attached.

143

The Green's functions G(L,L), G(L-1,L-l) and G(L-1,L) are
obtained through the recursion relationship described above.
5.3.1 Quantum percolation threshold:

The threshold for classical bond percolation is known to
be pc=0.5 for a two dimensional square lattice. The quantum
threshold pg is defined as the p value below which all
quantum states are localized. Obviously due to electronic
back scattering, it is expected that pg>pc. The scaling
theory of localization predicts that any finite degree of
disorder will result in all electronic states being
localized, i.e. pg=1.

It is noticed that if we apply the expression (5.55) to
calculate conductance, we always obtain finite values even
for a perfect lattice. Our calculation shows that g(LL)=Lq
the system width. Thus, pg is determined from calculating
the conductance as a function of system size L at a fixed
disorder quantified by p. We performed the calculation of
g(p,L) as a function of L at fixed p values(see Figure 5.2).
The results show that if p is very close to 1, g(p,L)
increases with L for small system size and reaches a maximum
at L=LM, and starts to decrease monotonically with L for
L>LM. When p is not very close to 1(p=0.95 or less), g(p,L)
is simply monotonically decreasing with the system size that
we chose. In other words, Lu decreases as p increases. The

point is that even for a very small disorder( p=0.98),

g(red. unit)

103(3)

40

3.7

3.4

144

 

F T I I U U V l I U I Y I V

(a): p=0.98

 

l I l l

L]

1

 

TUITU'UIUVUVITU'TIUTUVIUU

 

d
d.

  

LlLllllll

 

Llllljlllllll

 

 

Elam: The reduced conductance g as a function'of system
size. 8 different sizes are considered. p is fixed at 0.98. Each
point is an average over 200 configurations. The conductance g
increaseswithLfirstandreechesanmtimalvalueatL=280,then
startstodeerease.lteon'eepondstothecrossoverfroma
"scattering" region to the weak localization region.

g(red. unit)

103(3)

 

145

 

16

10

q
«4
q
d

—(
d

V r I l t I r

(a): p-O.97

d
_.(
q

 

3.1

rIUIVIIVlTYfIrU'TUIUfT‘VI'IUtj‘I'IUTVl'U'rrIUT‘

 

db

1.
2b....“o.

J)-
..
..
..e
3'
..
.(

(b): p=0.97

JlLllLlLllllLllJllllllllllLllLllllLlllllllll

 

 

2.7
0

100 200 300

W A similar plot to Figure 5.2 except that the calculation
is done for site percolation. The p value is fixed at 0.97 and again
each point is obtained from taking an average over 200 samples. g
increases with 1.. first but it reaches a maximum at 1.8120 and
starts to decrease.

 

146

g(p,L) will eventually decreases with L monotonically.
Eigg;g_§‘z,shows the numerical results of g(p,L) of a bond
percolating lattice. Each point in the plot is an average
over 200 configurations. A similar behavior is seen in
zigg;g_§;1 which is the result from calculations carried out
on a site percolating square lattice. As will be shown in
the next subsection, g(p,L) decreases exponentially with L
when the system size is sufficiently large, for finite p.
Therefore, it is clear that the quantum threshold for a 2-
dimensional square lattice is 1; any finite disorder(site or
bond) will result in the localization of electronic states.
5.3.2 Scaling of conductance:

The heart of the scaling theory of localization is the
one-parameter assumption that the scaling function depends
on 9 only. In this subsection, our objective is to perform a
detailed calculation of g(p,L) and the B function and to
analyze their behaviors in several distinct regions.

£1gu;g_§L1 shows a log-linear plot of g(p,L) as a
function of p. g(p,L) monotonically decreases with
decreasing p, it is not difficult to divide the figure into
three regions. In region (I): In g(p,L) decreases rapidly, a
very small decrement in p results in a big change in 9. We
call this the scattering region because when (1-p) is small
the defects are very sparse and the mean free path exceeds

the system size L. In region (II): g(p,L) starts to degrade

147

 

 

 

 

 

 

 

 

 

 

 

I : I l : I 1 I :
012 0!:
(b):L-32 ‘
100 _
' I
m-a L—
10'4 1
0.0

 

MAlog-lineerplotofgasafimctionofpatfixed
system size. Three distinct regions are classified: region (I) is very
narrowinwhichgdeeaysveryfastiterossesovertoregion(ll)
at p=0.05. Regional) is the weak localization region in which g
decays relatively slowly and g decreases as 1. increases. As more
bonds are removed, g enters the strong loealization region, region
(III) in which g decays much fasterthan in region (11).

148

slowly and it is called the weak localization region.
Finally in region (III): g(p,L) decreases very rapidly once
again. This corresponds to the strong localization region.
Our analysis is based upon the calculations of g(p,L)(or

ln g(p,L)) in those regions. Let pc1 and pcz be the points
separating those regions. In general, PCI and pcz depend on

system size. A more precise definition of pcl and Pcz is as

follow.
For a particular L, pcl is defined as
@(L’pcl)=0 (5.61)
éL

When the impurity concentration is very small, the mean free
path of electron scattering is much larger than the system
size L. It has been shown that a single defect will cause a
decrement in conductance by order of unity in that limit. By
noting that the conductance g is proportional to the system
size L when p=1, we have the following scaling form:-

g(P.L)~L-a(1-p)L2 (5.62)
with parameter a being of order of one. Therefore Pbl can be
roughly obtained as a function of L:

l-pc,~l/L (5.63)
as L-NO, pcl-al. In the region p>pc1, g increases as L
increases, so B is positive, which contradicts the scaling
theory of localization, which predicts non-positivity of B
in 2-dimensional disorder systems. It should be pointed here

that the scattering region is not physically reasonable as

149

it fails to yield an infinite conductance for a perfect
system.

In the weak localization region, g is monotonically
decreasing with L, and the B function is negative. We varied
the p values and calculated the conductance g as a function
of system size. An attempt has been made to fit the data to
the equation g(L)=g0-cln(L/L0) and the data fits the curve
reasonably well as shown in Eiggze 5,5. It implies that the
B-function scales as

fl=—a/g (5.64)
Where the parameter a is equal to the slope of the line in
Eigg;g_§‘§. The result is perfectly consistent with
perturbation theory.

In region (III), g(L) crosses over to a strong
localization region. Care must be taken when calculating the
conductance as g becomes extremely small when p approaches
pc, so that the n factor in the Green's function(Equation
(5.49)) starts to play a significant role in the numerical
calculations. Generally, small 9 requires small n. Eiggzg
5‘6 presents the results on conductance for relatively large
disorder. The logarithmic function of g(p,L) was plotted
against system size L, showing the linear dependence between

them. The slope of the line is dependent on disorder alone.

2.0

1.8

1.6

1.4

1.2

1.0

0.8

150

 

[Tllrrllll1rrlrllllrlrr

I'll

 

-1
—

r

l I IIIITII 1

LA

iliiiilrntlirirltLulir

 

—_
l...
b
h
l.
P
i—
-

 

50 100 200 500
103(1)

p...
O
N .
O

W Linear-log plot of g as a function of system size L,
showingtbebehaviorofgintheweakloealization region. Inthe
graphlhavetakenpao.075andthesysternsizesarechosenz
L816, 32, 64, 128, 256. The data for L=256 is from 100
configurationsandtherestofdataarefmmaZOOsample
average. It shows nearly a perfect linear relationship.

151

 

O
__
..
..
..
..

    
  
   

 

IIIIIIWIIIIITIIjI

 

 

: p=0.60.r=4.0 I
1:. -10 ‘L" —.
‘67) r ..
.9. C :
-15 :- =0.60.r=5.0 —.
- .4
Pl 1 1 1 1 l4 1 1 1 l 1 1 1 J l 1 1 1 1 J 1 1-

50 100 150 200 250

size L

Mi A log-linear plot of g as a function of L. It
corresponding to the strong localization region. The calculation
wasdoneforthesitepercolationeasewherep istakenas 0.60. r
stands for the energy level for an absent site while r=0 for any
present site. The upper curve has a little deviation from linear
dependencewhilethelowercuryeshowsanice straight lineas
predicted from strong localization.

 

152

5.3.3 Conductance distribution:

The distribution of conductance has been investigated
first by B. Sharpiro(Sharpiro, 1986,1987) and coworkers
based on the concept that two quantum resistors p, and ,02 in
series have the total resistance:

p=p. +102 +2lp.p2(1 +p. )(1 +7001"2 c089= "(13.49560 (5 . 55)
where 9 is the phase angle between )0, and p2. Suppose the
distributions of p, and p2 are known and denoted by P(p) and
P(p) , the distribution of p under condition p=u(p,,p2,€) is

given by:
l
P(p)=;Jdp.dp.d6P.(p.)I:(p.)6<p—u) (5...)

Assume that pI and p2 have linear scales AL and L
respectively with AL being small and thus small p1, one
expands the above expression up to linear terms of pi,
yielding:

M] (5...)

PL+AL(p) = PL(p)+(/)l>[§)(p+p2) 5p

where (,) is the average value of p1 and has been assumed to

be linearly dependent on AL as it is small. Let a=(p,)/AL,

one obtains a diffusion like equation:

OPLU’L j 2 at100)
d -a{ap(p+p) 5p] <5-68)

PL(p) can be easily solved for two limits of pw+0 and p—wn.
The former limit implies that the p2 term can be neglected

and thus the solution is written as:

P.(p)=iexp(-p/a£) (5.59)

153

The average of p yields <p>=aL, recovering Ohm's law. When p
>>1, p2+p5p2, (5.68) turns into a diffusion equation and the

solution reads:

1 .
PL(p)=mexp[—(lnp—aL)‘/4aL] (5.70)

In this case, the average of lnp depends linearly on L. Both
distributions support the single-parameter scaling theory
assumption.

The conductance distribution in two dimensions is
difficult to find because of the traverse scattering. Many
techniques such as the "fan transformation"(P.W.Anderson,
1980) and a Kadanoff type scaling transformation are
employed to generalize the one dimensional result. Here we
directly calculated the distribution of conductance and
found the distribution is log-normal when g is very small as
expected. 2199;; 5,7 shows a typical distribution for a
strongly disordered system and the curve fits a gaussian
distribution well. But for medium disorder, the distribution
is anomalous and of bimodal structure. Eiguze 5.8 is the
distribution for a 16x16 square lattice with p=0.75. The
statistics was done for some 4000 samples. Clearly, the
graph shows that P(g) has two maxima: one is located near
g=1 and the other is close to g=0. When disorder is small,
the peak near g=0 is not seen. The peak near g=1 decreases

as disorder increases and will eventually disappear,

154

 

 

 

 

 

l I
e
400 — °
e e.
e
e
e
e
200 i—
e
e
e
e
. e
e
e e
W. l o’fJ
-40 -20 O

M Distribution of g: strong localization limit. The
statisticswasfor4000configurations.Lisfixedat32,andpis
fixedatO.40.Thex-axisislog(g)andycoordinaterepresentsthe
nurnber.Thepealtisloeatednesrlog(g)-15.'I'heoverallshape
eanbefittedcloselytoagaussianeurve.

 

155

 

(D
O
O
.1
.1
J
—(
..l
..._l
1
.J
—(
1

 
   

400 ‘

300

 

200

LlllllllLlLLllIlJllIlLlllllLJ

Ililrlllrlrriill1ti

 

 

 

00
ol-
)-
L.
P—
on
)-
l-
L.
3"—
oli-
2‘1.
01
.N
o
9"
on

Wbisn'ibunonofg:mediumdegieeofdisolder.fleredie
sample number is 4000 again. system size L-16. p-0.25. I have
takenafinitehoppingconstantx-OJ insteedofOforanabsent
bondtoavoid floatingnumberunderflow. ltisclearthsttwo
mkmshowup:cneisnesr0,theotherisneerg-I.O.The
bimodalfestureleedstoadifficultyintakingavemges.lflog(g)
isavengeithecontributionfi'omgneerOdominatesflfgis
averaged.thecontributionfromg-l dominates.

156

merging into the strong localization case. So far, we do not

have physical explanation for this bi-modal structure.

u s .

It is interesting to study the scaling behavior of
conductance of a lattice near its classical percolation
threshold pc, where the backbone of the lattice structure is
fractal. It has been discovered that diffusion on a fractal
structure exhibits abnormal behavior and the diffusion
constant D algebraically decays with distance r as discussed
in §gg§ign_551‘z. About half a decade ago, several
authors(Shapir and Aharony,1982;Levy and Souillard,1987;
Lambert and Hughes,1991) studied quantum diffusion on a
fractal structure using scaling arguments and introduced the
concept of super-localization in contrast to the
localization behavior discussed earlier. Super-localization
claims that the electronic wavefunction on a fractal
structure decays faster than exponential against distance
and typically satisfies the following scaling form:

MVP M0)CXPl-(r/ 6)“) (5-71)
with a>1. This implies that the conductance g of a fractal
structure has a similar scaling with its size L. Because the
conductance is directly related to the electronic
wavefunctions, it is very desirable to carry out a

systematic numerical calculation of conductance to verify

157

this relationship. It is noticed that near the percolation
point, the conductance is extremely small due to fewer
connected points, so in order to obtain reliable calculation
results, the hopping elements in the Hamiltonian are set to
be either 1 or a small finite number. This procedure will
certainly effect the consequences of calculations but it
does not effect the scaling behavior as no structural change
has been made. The two dimensional lattice sizes are
selected to be L=16, 32, 48, 64, 80, 96,112, 128, 144,160.
The conductance log(g) as a function of L is plotted in
£1gg;g_552, showing a linear dependence as predicted from
strong localization. It should be pointed out that because
of the introduction of a small hopping constant, the result
shows a normal localization behavior. More detailed
calculation needs to take the small hopping constant as
small as possible. Unfortunately, our method cannot handle
the situation well due to computer rounding error. One
appealing method is to utilize the transfer matrix approach
discussed in ggggign_552. It remains a prospective project
in the near future.

To conclude, we have calculated the conductance of a 2d
percolation lattice as a function of system size and p
value. Our calculation results strongly support the scaling

theory of localization.

 

conductance (log(g))

158

 

 

 

 

0.0 _ . ,
i 2
-2.5 :- —’
L I
i- a
-5.0 :- -—.
E 2
-7.5 _— —~
-1o.o — —"-
C I
-12.5 :- -:
-15.o ' . . ‘
O

 

Ming-linearplotofgasfimctionofLatpercolation
point p-0.5. To avoid singularities, the hopping constant for an
absentbondhasbeentakentobeO.l.Eachpointinthefigureis
an average of 400 configurations. We wish to see a super-
loeslization behavior equation (5.71). Instead, it shows clearly a
lineerdependence,aspredictedfrornstronglocalization.

159

21551_ggm;;55; In the final section, I will provide a
limited discussions on three problems related to the work in
the thesis.

(a) Modification to the ICU model:

In the interrupted coalescence model(ICM), two
unrealistic assumptions have been made: (i) it is assumed
the surface tension of the droplets is isotropic so that the
equilibrium shape is a spherical cap. (ii) no dynamic
process has applied to the two large touched droplets whose
coalescing process has presumably been interrupted. The
touched area is unphysical as there exist singular points.
Both of those assumptions contribute to the deviation of the
simulation figures from the experimental results(fiigggg
551). Careful examination of £1gg;g_551,shows that the
droplet shape is not spherical, instead several well
developed facets are seen even in the early stage where
complete coalescence takes place. The complete coalescence
time in an anisotropic model has a stronger dependence on
system size than that in an isotropic model, which implies
the size effect may play an important role in the
interruption mechanism. The touched area of two large
droplets has to be smoothened out. A rather rigorous method
of doing this is to introduce the surface diffusion process
as discussed in figggign_252. However it becomes impractical

as the number of those droplets is large. The calculation in

160

Chapter Two shows that the healing process in the neck area
is so speedy that only a small fraction of the total
coalescence time is needed to reach the neck elimination
stage. Therefore, as an approximation, the smoothening
procedure may directly incorperate the neck elimination
stage. The simulated figure would then be much closer to the
SEM photograph.

(b) Dynamic phase in Au clusters:

A solid phase has different dynamic behavior from a
liquid phase and whether a cluster is in a solid phase or in
a liquid phase depends not only on the temperature but also
on the cluster size. Our simulation results show that a
solid cluster may exhibit distinct dynamic behaviors(£1gg;§
1511) at various sizes and temperatures. The physical basis
for the existence of those different dynamics is not very
clear. Our tentative explanation is that diffusion which
occurs in large solid clusters but is absent in small solid
clusters.

It has been pointed out that the strong existence of
surface diffusion in solid gold clusters is due to the many
body interaction which suggests a stronger coordination
dependence than pairwise interaction(Ercolessi etc.,1988).
An obvious choice is to test a Lennard-Jones(LJ) system and
observe if the system has a similar behavior to gold

clusters. Many MD calculations show that the system of LJ

161

interaction exhibits a solid-liquid phase transition. The
surface atoms in any solid cluster have less binding energy
than their bulk counterparts and therefore may also have
similar dynamic behavior as observed in gold clusters. A
detailed MD simulation on LJ systems will answer part of the
question.

(c) Super-localization problem:

Scaling arguments reveal that due to the fractal
dimensionality of the incipient infinite cluster of a
percolation structure near its threshold, the diffusion on
such a structure is anomalous. Both theoretical arguments
and numerical simulations showed that the phonon's amplitude
decays faster than exponential as stated in equation(5.71)
on fractal structures. It is very desirable to perform a
detailed numerical calculation of conductance of a system
with fractal dimension. As pointed out in Chapter 5, the
conductance near threshold is so small that a finite hopping
constant(ideally zero) for an absent bond has to be used to
obtain reliable calculation results. The introduction of
this small number may profoundly alter the nature of its
wavefunction, so one has to make the value as small as
possible. The transfer matrix method described in sgction
552 provides an alternate and perhaps more efficient way to
calculate conductance. Along with the calculation, two

problems have to be addressed: (i) the relationship between

 

162

the scattering matrix and conductance, i.e. the extension of
Landauer's formula; (ii) the conversion from the transfer
matrix to the scattering matrix.

As a closing statement, I quote the words from Einstein
"The nature is as complex as it needs to be, no more".
Although the nature is asymmetrical and dynamic( it is often
called beauty and mystery), life would be much easier if it
is understood in its least complex form, which may be the

ultimate goal for physicists.

APPENDIX A

TUVYTVTTIITTYITUTITITYTWTTYIITTIYTTTTTITrTTTTTTITITTTTTITTTITT

 

C1111111L1111111L1111111111L1111111L11111111111111111111111111L1
T

The program is created to study properties

C
C

of Au-cluster via constant temperature MD

CL111111111111111L111111111111111111

LLLLL

AALLJLLLLILLLILJlllLJLILA

 

TTTU VTTI IIWVIUTYTTYWIYIWWTTTIVTT' VjTTYTTTT‘rTTTTTTTITTIITTTTTIT

LllllllLlnLl1LLLLJALIJLLLLJLJLllLLl11lLllJLIJllLLLLllALLAILl11

(:TTTTTITTITIWY'WTTTIVY ITTYTWITITTITTITTTTTITTTW'IT1ITTTTIITTTT

 

Subroutines l.readin

 

.ww n9
mflfinnnnfl
I.nfOs&Wfl
.m..no..c..h. .e.
2345678

C
C
C
C
C
C
C
C

ILLILLALLAILILI1111ALILLILLLJllLllLIJLJsllliLllLLLLLLJILLLLLLLL

(:TTTYITTIIITTTUTTTTIYYTTTYIFTTTWTTYTTYTITTTTTTTITTTWTTTTIIIITT

 

parameter (n=6l3)

,vx(800),vy(800),vz(800)
,ax(800),ay(800),az(800)
,bx(800),by(800),bz(800)
,cx(800),cy(800),cz(800)
,fx(800),fy(800),fz(800)

common/circles! rx(800),ry(800),rz(800)

&&&&&

common/circles/ s,sv,sa,sb,sc,sf

common/circles/ vexp(20000), vexq(20000)

common/circles/ nei(90000), lpoint(800)

integer free

p=10.15
q=4.l3

aa=0.ll8438

dt = 0.02
ne = le4

(M = 2-*q

Onewl)

open (unit=29, file='rg.out', status

100

164

open (unit=28, file='[yuxinhua.md]en.out', status='new')
nstep = 2e5

ivt = le2

tmp = 1.0e-3

rcut2 = 3.0

qh = 40.

free = 3

call readin ( n,dt )

do 100 ie = l,2“‘ne

vexp(ie) = exp((] .-real(ie-0.5)/real(ne))"‘p)
vexq(ie) = exp(( 1 .-real(ie-0.S)/real(ne))*qq)
confinue

nfree = n*free

do itmp = 1,6
tmp = 5.0e-3 - (4.0e-4) "' float(itmp)

ten = 0.
Mm=0.

do 1000 istep = l,nstep

int = mod(istep-l,20)

if ( int.eq.0) call listn ( n, rcut2)

call predic (n, dt)

call force ( n, (it, aa, p, q, ne)

call kinet ( n, en_k )

sf= (2.0 " en_k- real (nfree + 1) * tmp ) / (s‘qh)
call correc (n, dt)

rxc = 0.0

ryc = 0.0
rzc = 0.0

 

1001

1010

1 000

165

do 1001 i= l,n
rxc =rxc+rx(i)

0'0 = We + r5'0)
rzc = rzc + rz(i)
confinue

rxc = rxc/real(n)

ryc = ryclreal(n)

rzc = rzc/real(n)

rg = 0.0

do 1010i = 1,n

r8 = r8 + (rX(i) - rXC)"2 + (Mi) - WC)"2 + (rZ(i)-r26)“2
confinue

rg = sqrt(rg/real(n))

arrargi-rg
if (isteple. 1e5) goto 1000

md=mod(istep,ivt)

if (md.ne.0) goto 1000
itime = (istep-1e5)fivt
call energy ( n,en,temp )
write ( 28,‘I )tmp,temp,en
ten=ten+en
atm=atm+temp
arg=arglreal(ivt)

write ( 29f) itime,arg
arg=0.

continue

ten=ten/ 100.

atm=atm/100.
write ( 28,‘ )tmp,ten,atm

 

166

end do
call writin ( n )

close ( unit=28 )
c close ( unit=29 )

stop
end

C +++-+-+-H—+-+-+—+-+—+-+-H- end of main program +++++-l-+++-+++—H—+—++-l—l—i—l-—+—
CH4+H++++H—FF+H—+++4—HH+F++FH+FH+++++H+++4—H4+H+H+++H++
C begin subroutine readin

CWWWWHWWH

subroutine readin (n,dt)

common/circles/ rx(800),ry(800),rz(800),

& vx(800),vy(800),vz(800),
& ax(800),ay(800),az(800),
& bx(800),by(800),b2(800),
& cx(800),cy(800),cz(800),
& fit(800),fy(800),fz(800)

common/circles/ s,sv,sa,sb,sc,sf

open ( unit=21, file='[yuxinhua.md]xyz.in', status='old')
open ( unit=22, file='[yuxinhua.md]vcc.in', status='old')
open ( unit=23, filfi'[yuxinhua.md]acc.in', status='old')
open ( unit=24, file='[yuxinhua.md]bcc.in', status='old')
open ( unit=25, file='[yuxinhua.md]ccc.in', status='old')
open ( unit=26, file=’[yuxinhua.md]scc.in', status='old')

do 100 i=l,n
read (21,"' ) a,b,c
rx(iFa
ry(i)=b
rz(i)=c

100 confinue

do 200 i=l,n
read (22,*) ab,c

 

200

300

400

500

C+-H—+—H+l-+-+-l- end ofreadin +++++§ H+1H *

167

vx(i)=a
vy(i)=b
vz(i)==c

confinue

do 300 i=1,n
read (23,"') a,b,c
ax(i)=a

ay(i)=b

512(ti

continue

do 400 i=1,n
read (24,*) a,b,c
bx(iFa

by(i)=b

bz(i)=c

confinue

do 500 i=1,n
read (253') a,b,c
cx(i)=a

cy(i)=b

cz(i)=c

confinue
read (26,"') s,sv,sa,sb,sc,sf

close ( unit=21 )
close ( unit=22 )
close ( unit=23 )
close ( unit=24 )
close ( unit=25 )
close ( unit=26 )

return
end

C+—+—+-+-+—+-+++—+—++-+-+ beginning of subroutine force +—H—-l—l-l-+-l—l—l—H—l—+-H++'l—+

 

eeeee

168

subroutine force ( n,dt,aa,p,q,ne )

common/circles/ rx(800),ry(800),rz(800),
vx(800),vy(800),vz(800),
ax(800),ay(800),az(800),
bx(800),by(800),bz(800),
cx(800),cy(800),cz(800),
fx(800),fy(800),fz(800)

common/circles/ s,sv,sa,sb,sc,sf
common/circles/ vexp(20000),vexq(20000)

common/circles/ nei(90000),lpoint(800)

qq=2-“q
rc2=2.25

do 200 i=1,n

rxi=rx(i)
ryi=ry(i)
rzi=rz(i)
fiti=0.
fyi=0.
fzi=0.

nlow=lpoint(i)
nup=lpoint(i+l)-l

rxp=0.
do 10 list=nlow,nup

j=nei(list)

rxij=rxi-rx(i)

ryijsyi-ryfi)

rzij=rzi~rz(i)
fi12=rxij*ndl+ryii’ryij+rzii*rzij
if(rij2.gt.rc2) goto 10

rij=sqfl(n'J°2)

 

10

180

169

ie=rij‘ne

rxp=rxp+vexq(ie)
continue

nor=sqrt(rxp)

do 190 list=nlow,nup

sxp=0.

j=nei(list)

do 1 80 ljst=|point(j),lpoint(j+1)-1
k=nei(ljst)

rxjk=rx(j)-rx(k)
rxik=w(i)-ry(k)
rzjk=rz(j)-rz(k)
Ijk2=rxjk"2+ryjk"2+rzjk"2
if(tjk2.gt.rc2) goto 180
Iik=sqn(tik2)

ie=ne"rjk

sxp=sxp+vexq(ie)

confinue

sxx=sqrt(sxp)
rxiiji-rx(nei(list))
ryii=ryi-ry(nei(list))
rzij=rzi-rz(nei(list))
rij2=rxij"2+ryij"2+rzij"2

if(rij2.gt.rc2) goto 190

rij=sqn(rb'2)
ie=ne*rij

 

170

fl=aa*p"vexp(ie)
12=q"vexq(ie)/r>o:
f3=q‘vexq(ie)/sxx

f12=fl -0.5"(f2+f3)

fid=fid+(rxij/rij)"fl 2

fyi=fyi+(ryil/rij)‘f12

fzi=fzi+(rzij/rij)"'f12
190 confinue

fx(i)=fxi
60)“in
fz(i)=fzi

200 confinue

return
end

C-+-++++-+-+-+-H-+-++H- end of subroutine force +H—i-+++-+-+—l—+—+-+-+-+—+—+—+—+-++-+—+-l—l-+

c++-+-+—+—+-+—+—++-+-+-+—+—++ beginning of subroutine kinet ++++++++++++++++++
subroutine kinet (n,en_k)

common/circles/ rx(800),ry(800),rz(800),
vx(800),vy(800),vz(800),
ax(800),ay(800),az(800),
bx(800),by(800),bz(800),
cx(800),cy(800),cz(800),
fit(800),fy(800),fz(800)

858°2°R°R°

common/circles/ s,sv,sa,sb,sc,sf

en-k=0.

do 100 i=1,n

en_k=en_k+vx(i)"2+vy(i)**2+vz(i)**2
100 confinue

en_k=(0.5"en_k)*(s*s)

 

171

return
end

C +-+-+-+-+-+-i-+-+—l-i-++++ end of subroutine kinet ++++H-l-i-+++-++++-l-H—h+++-+—i—+—+

C+—+-+—+-+—+—+++++-l-+-+-+ beginning as predictor +—+—++l—++H—+-i—+++-+—+-+-+—+-+-+++-l-l--+—+-
subroutine predic ( n,dt)

common/circles/ rx(800),ry(800),rz(800),

& vx(800),vy(800),vz(800),
& ax(800),ay(800),az(800),
& bx(800),by(800),bz(800),
& cx(800),cy(800),cz(800),
& fx(800),iy(800),fz(800)

common/circles/ s,sv,sa,sb,sc,sf

cl=dt

c2=cl *dt/2.0
c3=c2*dt/3.0
c4=c3"'dt/4.0

do 100 i=1,n

rx(i)=rx(i)+cl *vx(i)+c2*ax(i)+c3 *bx(i)+c4*cx(i)
l'>'(i)=ry(i)+cl ‘Vy(i)+02*ay(i)+c3 *by(i)+c4*cy(i)
rz(i)=rz(i)+cl *vz(i)+c2*az(i)+c3 ’bz(i)+c4*cz(i)

vx(i)=vx(i)+cl *ax(i)+c2*bx(i)+c3 I"cx(i)
W(i)=VY(i)+cl ‘aY(i)+02‘by(i)+c3 *CYG)
vz(i)=vz(i)+c1 ’az(i)+c2"‘bz(i)+c3 *cz(i)

ax(i)=ax(i)+cl ‘bx(i)+c2"'cx(i)
ay(i}=ay(i)+01 *by(i)+02"cy(i)
az(i)=az(i)+c l ‘bz(i)+c2"‘cz(i)

bx(i)=bx(i)+c1 1’cx(i)
b)’(i)=b)'(i)+cl ‘CY(i)
bz(i)=bz(i)+cl *cz(i)

 

172

100 confinue

s = s+c1 ‘sv+c2"sa+c3*sb+c4*sc
sv=sv+c1 " sa+c2* sb+c3 " sc
sa=sa+cl "' sb+c2“' sc

sb=sb+cl "sc

return
end

C++—+—+-++-i—+—+—+—+-l—+-+-+—+—+—+—+—l— end of subroutine predictor ++-+—1—+—+—+—+—+-+++-+++++-+—++

C++—+-++-+—+-+-+—+—+—+-++++—+-+—l- beginning of subroutine corrector +-+-+-+-i-+++++-H-++

subroutine correc (n,dt)

common/circles/ rx(800),ry(800),rz(800),
vx(800),vy(800),vz(800),
ax(800),ay(800),az(800),
bx(800),by(800),bz(800),
cx(800),cy(800),cz(800),
fir(800),fy(800),fz(800)

asses

common/circles/ s,sv,sa,sb,sc,sf

parameter ( g0=l9.0/90.0, g1=3.0/4.0,g3=1.0/2.0,g4=1.0/ 12.0 )
cl=dt

c2=cl*dt/2.0

c3=c2"'dt/3.0

c4=c3"‘dt/4.0

cr=g0*c2

cv=gl "c2/cl

cb=g3"c2/c3
=g4‘c2/c4

dmmem
'axi=fx(i)/(s*s)-2.o*sv*vx(i)/s

ayi=fy(i)/(s‘s)-2.0*sv‘vy(i)/s
azi=fz(i)/(s*s)-2.0*sv*vz(i)/s

 

173

corx=axi-ax(i)

cory=ayi-ay(i)
corz=azi-az(i)

rx(i)=rx(i)+cr*corx
ry(i)=ry(i)+cr‘cory
rz(i)=rz(i)+cr*corz

vx(i)=vx(i)+cv*corx

vy(i)=vy(i)+CV*cory
vz(i)=vz(i)+cv*corz

ax(i)=axi
ay(i)=ayi
az(i)=azi
bx(i)=bx(i)+cb"corx

by(i)=by(i)+cb“‘cory
bz(i)=bz(i)+cb*corz

cx(i)=cx(i)+cc*corx
cy(i)=cy(i)+cc*cory
cz(i)=cz(i)+cc*corz

100 confinue
cors=sf-sa
5 =5 +cr‘cors
sv=sv+cv*cors
=sf
sb=sb+cb"cors
sc=sc+cc"‘cors

return
end

C W end of subroutine correc -+-++-+—+—H—+—+—i—+—+—+-+-i—+—H-+—+--+—l—+

C++++-+—+++—+—+—++ beginning of subroutine listn +-++-l—+—+-+-++-l—+-+-+-+-+-+-i-+-+—+-i—+-

subroutine listn (n,rcut2)

 

sesss

199
100

174

common/circles/ rx(800),ry(800),rz(800),
vx(800),vy(800),vz(800),
ax(800),ay(800),az(800),
bx(800),by(800),bz(800),
cx(800),cy(800),cz(800),
fit(800),fy(800),fz(800)

common/circles/ s,sv,sa,sb,sc,sf
common/circles/ vexp(20000),vexq(20000)
common/circles] nei(90000),lpoint(800)
integer list
list=0
do 100 i=1,n
rxi=rx(i)
wi=ry(i)
'=rz(i)
lpoint(i)=list+1
do 199 i=1,n
iflieqj) goto 199
rxij=rxi-rx(j)
ryii=ryi-ry(i)
rzij=rzi-rz(j)
rij2=ntij”2+ryij“2+rzij**2

if(rij2.le.rcut2) then

list=list+ 1
nei(list)=j
endif
confinue

continue
lpoint(n+ 1)=list+ l

 

175

return
end

C++++-l-H—+—+—+-+++-+—+-i- end of subroutine listn ++-+—+—+—+—+—+—+—+—+—1-—++++++++++++—+—+—+

c+-+—+-+-l-+—+-++~+—+—+-+—+-++ beginning of subroutine writin +++-+++-+-++++++++++++++

subroutine writin (n)

common/circles/ rx(800),ry(800),rz(800),
“(800).W(300).VZ(800).
ax(800),ay(800),az(800),
bx(800),by(800),bz(800),
cx(800),cy(800),cz(800),
fx(800),fy(800),fz(800)

sssss

common/circles/ s,sv,sa,sb,sc,sf

open (unit=2 l ,file='[yuxinhua.md]xyz.in',status='old')
open (unit=22,file='[yuxinhua.md]vcc.in',status=’old‘)
open (unit=23,file='[yuxinhua.md]acc.in',status='old')
open (unit=24,file='[yuxinhua.md]bcc.in’,status='old')
Open (unit=25,file='[yuxinhua.md]ccc.in',status='old')
open (unit=26,file='[yuxinhua.md]scc.in',status='old')

do 100 i=1,n

write (213') rx(i),ry(i),rz(i)

write (223) VX(i).Vy(i).VZ(i)
write (211,“) ax(i),ay(i),az(i)
write (24,"') bx(i),by(i),bz(i)
write (253') cx(i),cy(i),cz(i)

100 confinue
write (26,"‘) s,sv,sa,sb,sc,sf

close (unit=21)
close (unit=22)
close (unit=23)
close (unit=24)
close (unit=25)

 

176

close (unit=26)

return
end

C ++—+—++—+—+++-+—+—+++++ end of subroutine writin +-+-+-+++-++-+—+—i—+—l—-+—H—+-+-i—+—+++-+—+-

c++++++-H—+—++++-+—+—++ beginnig of subroutine energy +4++-+—-+—+—+—+—+-+—++—+++—+—+—+++

subroutine energy (n,en,temp)

common/circles/ rx(800),ry(800),rz(800),

& vx(800),vy(800),vz(800),
& ax(800),ay(800),az(800),
& bx(800),by(800),bz(800),
& cx(800),cy(800),cz(800),
& fx(800),fy(800),fz(800)

common/circles/ s,sv,sa,sb,sc,sf

common/circles/ vexp(20000),vexq(20000)
common/circles/ nei(90000),lpoint(800)

temp = 0.0

do 100 i = l,n

temp = temp + (vx(i)**2 + vy(i)"2 + vz(i)“2)
100 confinue

temp=0.5*(s"s)‘temp

p= 10.15
q=4.13
aa=0.118438

qq = 2-‘q

vn = 0.0

do 200 i=1,n
vn1=0.
vn2=0.
rxi=rx(i)
ryi=ry(i)
rzi=rz(i)

 

177

do 190 j=1,n

if (j.eq.i) goto 190
ncij=nd-rX(i)

ryij=ryi-ry(i)

rzij=rzi-r2(i)
rij2=rxij"2+ryij"2+rzij"2
rij=sqrt(fi1'2)

vn1=vnl+exp(-qq*(rij-l))
vn2=vn2+exp(-p*(rij- 1))

190 confinue

vn1=sqrt(vn1)

vn=vn+0. 5 ‘(aa‘an-vn 1)
200 confinue

en=temp+vn

en=enlreal(n)
temp=2. "temp/rea1(3“n)

return
end
C-l—+-+-++++—i-+-+-+-l—+—l-++-+-l—+—l- end of energy +1-H—H—H—i—i—l—i-++-l—l-+-i—+—+++++++++

 

W
1. Abrahams, E., Anderson, P. W., Licciardello, D. C. and Ramakrishnan, T. V.,

Physical Review Letters 4;, 673 (1979).
2. Abraham, F. F., Advances in Physics 3, 1 (1986).

w

Ajayan, P. M. and Marks, L. 0., Physical Review Letters @, 585 (1988).
Ajayan, P. M. and Marks, L. D., Physical Review Letters 51, 279 (1989).
Andersen, H. C., J. Chem. Phys. 12, 2384 (1980).

Anderson, P. W., Physical Review 113, 1492 (1958).

Anderson, P. W., Physical Review _1_2_4, 41 (1961).

Anderson, P. W., Physical Review B 2_3_, 4828 (1981).

§OP°>’.O‘S":‘>

Andersen, P. W., Thouless, D. J., Abrahams, E. and Fisher, D. S., Physical Review B

22, 3519 (1980).

10. Avron, J. E., Taylor, J. E. and Zia, R. K. P., Journal of Statistical Physics _3_3, 493
(1983)

11. Azbel, M. Y., .1. Phys. C 1;, L797 (1980).

12. Azbel, M. Y., Physics Letters 15; 410 (1980).

13. Azbel, M. Y., Physical review Letters 51, 1787 (1991).

14. Bergmann, 6., Physics Reports 191, l (1984).

15. Berry, R.S., in Few-body systems and Multiparticle Dynamics edited by D. Micha
(1987)

16. Bonzel, H. P., Preuss, E. and Steffen, 3., Appl. Phys. A 55, 1 (1984).

17. Bonzel, H. P., Preuss, E. and Steffen, 8., Surface Science 115, 20 (1984).

18. Borel, J.P, Surf Sci. L95 , l (1981)

19. Brailsford, A. D. and Gjostein, N. A., Journal of Applied Physics fl, 2390 (1975).

20. Briant, C. L. and Burton, J. J., J. Chem. Phys. 5;, 2045 (1975)

179

21. Briscoe, B. J. and Galvin, K. P., J. Phys. D _25, 1265 (1990).

22. Brooks, A. and Aharony, A., Europhys.Lett. fl, 1355 (1987).

23. Bufl‘at, Ph., Borel, J.P., Phys. RevA 1; 2287 (1976)

24. Bulgac, A. and Kusnezov, D., Physical Review Letters §_8_, 1335 ( 1992).

25. Buttiker, M., Physical Review Letters 51, 1761 (1986).

26. Car, R. and Parrinello, M., Physical Review Letters 22, 2471 (1985).

27. Castellani, C., Kotliar, G. and Lee, P. A., Physical Review Letters 2, 323 (1987).

28. Chui, S. T. and Weeks, 1. D., Physical Review B E, 4987 (1976).

29. Cleveland, C. L., J. Chem. Phys. 82, 4987 (1988).

30. Clogston, A. M., Physical Review 125, 439 (1962).

31. Cohen, A., Roth, Y. and Shapiro, 8., Physical Review B 3_8_, 12125 (1988-1).

32. Cowley, R. A., in ”Equilibrium Structure and Properties of Surfaces and Interfaces"
edited by Gonis and Stocks, Plenum, New York(1992).

33. Derrida, B., Godreche, C. and Yekutieli, 1., Europhys. Lett. 12, 385 (1990).

34. Dubson, M. A. and Jeffers, 6., Phys. Rev. B (in Press)(1994)

35. Dundurs, J ., Marks, L. D. and Ajayan, P. M., Philosophical Magazine A _51, 605
(1988)

36. Duxbury, P. M., Dubson, M., Yu, X. and Jeffers, G., Europhys. Lett. 2_6_, 601 (1994).

37. Economou, E. N. and Soukoulis, C. M., Physical Review Letters 4_6_, 618 (1981).

38. Ercolessi, F ., Andreoni, W. and Tosatti, 13., Phys. Rev. Lett. Q, 911 (1991).

39. Ercolessi, F., Parrinello, M. and Tosatti, 13., Phil. Mag. A g, 213 (1988).

40. Family, F. and Meakin, P., Phys. Rev. Lett. 51 , 428(1988).

4]. Family, F. and Meakin P., Phys. Rev. A Q, 3836 (1989).

42. Fisher, D. S. and Lee, P. A., Physical Review B 2. 6851 (1981).

180

43. Fisher, M.E., in Statistical Mechanics of Membrants and Surfaces edited by D.
Nelson, T. Piran and S. Weinberg(l988)

44. Friedel, J ., in The Physics of Metals edited by Ziman (1969).

45. Garzon, I. L. and Jellinek, 1, Z Phys. D 2_Q, 235 (1991).

46. Grilli, M. and Tosatti, 5., Physical Review Letters 52, 2889 (1989).

47. Gruber, E. E. and Mullins, W. W., J. Phys. 25, 875 (1967).

48. Gupta, R. P., Phys. Rev. B Q, 6265 (1981).

49. Harris, A. 8., Physical Review B 2_9_, 2519 (1984).

50. Herring, 0, Phys. Rev., 5; , 87 (1951).

51. Herring, C. in ”Structure and Properties ofSurfaces", edited by Gomer, R. and Smith,
C.S. (University of Chicago Press, 1952).

52. Hoover, W. G., Physical Review A .31, 1695 (1985).

53. lijima, S. and Ichihashi, T., Physical Review Letters 55, 616 (1986).

54. Jayaprakash, C., Saam, W. F. and Teitel, 8., Physical Review Letters 5_0, 2017 (1983).

55. Jeffers, G., Dubson, MA. and Duxbury, P. M., J. Appl. Phys. E 5016(1994).

56. Jellinek, 1., Beck, T. L. and Berry, R.S., J. Chem. Phys. E, 2783(1986).

57. Jellinek, J. and Berry, R. S., Physical ReviewA fl, 2816 (1989).

58. Jellinek, J. and Garzon, I. L., Z. Phys. D 25, 239 (1991).

59. Jiang, Z. and Ebner, C., Physical Review B 4_0_, 316 ( 1989).

60. Joanny, J. F. and Gennes, P. G., J. Chem. Phys. 51, 552 (1984).

61. Kellett, B. J. and Lange, F. F ., J. Am. Ceram. Soc. 12, 725 (1989).

62. Kosterlitz, J. M., J. Phys. C 1, 1046 (1974).

63. Lambert, C. J. and Hughes, G. D., Physical Review Letters Q, 1074 (1991).

64. Lancon, F. and Villain, J ., Physical Review Letters Q, 293 (1990).

181

65. Lancon, F. and Villain, in "Kinetics of Ordering and Growth at Surfaces" edited by
Lagally, M.G., Plenum, New York(l990).

66. Landau, L. D. and Lifshitz, E. M., Iheory of Elasticity, 3rd edt. (1986)

67. Landau, L. D. and Lifshitz, E. M., Electrodynamics of Continuous Media(1984)

68. Lange, F. F. and Kellett, B. J., J. Am. Ceram. Soc. 12, 735 ( 1989).

69. Langreth, D. C. and Abrahams, 5., Physical Review B 25, 2978 (1981).

70. Lee, P. A. and Fisher, D. S., Physical Review Letters fl, 882 (1981).

71. Lee, P. A. and Ramakrishnan, T. V., Reviews of Modern Physics Q, 287 (1985).

72. Lee, P. A., Stone, A. D. and Fukuyama, H., Physical Review B 55, 1039 (1987-11).

73. Levy, Y. -E. and Souillard, B., Europhys. Lett. 51, 233 (1987).

74. Matsuoka, H., Hirokawa, T., Matsui, M. and Doyama, M., Physical Review Letters 52,
297 (1992).

75. Mazor, A., Srolovitz, D. J., Hagan, P. S. and Bukiet, B. G., Physical Review Letters
@, 424 (1988).

76. Meir, Y. and Wingree, N. S., Physical Review Letters Q, 2512 (1992).

77. Miller, K. T. and Lange, F. F ., Acta Metall 51, 1343 (1989).

78. Mullins, W. W., J. Appl. Phys. 25, 333 (1957).

79. Mullins, W. W., J. Appl. Phys. 3_0_, 77 (1959)

80. Nichols, F. A. and Mullins, W. W., Journal of Applied Physics 55, 1826 (1965).

81. Nose, S., J. Chem. Phys. 5_l_, 511 (1984).

82. Nozieres, P. J. Physique _45, 1605 (1987).

83. Ozdemir, M. and Zangwill, A., Physical Review B _4_2_, 5013 (1990-1).

84. Pandit, R., Schick, M. and Wortis, M., Physical Review B 55, 5112 (1982).

85. Pashley, D.W. Adv. Phys. E , 327(1965).

 

182

86. Pashley, D. W., Sowell, M. 1, Jacobs, M. H. and Law, T. J., Philos. Mag. l_9, 127
(1964).

87. Preuss, E., Freyer, N. and Bonzel, H. P., Appl. Phys. A 4_1 , 137 (1986).

88. Raghavan, R. and Mattis, D. 0, Physical Review B 25, 4791 (1981).

89. Robbins, M. O. and Joanny, J. F., Europhys. Lett. 5, 729 (1987).

90. Schommers, W. and von Blanckenhagen, P., "Structure and Dynamics of Surfaces"
(1987)

91. Selke, W. and Duxbury, P. M., Z. Phys. B _9_4_, 311-318 (1994).

 

92. Shapir, Y. and Aharony, A., Physical Review Letters Q, 486 (1982).

93. Shapiro, 8., Physical Review B 55, 4394 (1986).

94. Shapiro, B., Philosophical Magazine B 55, 1031 (1987).

95. Soukoulis, C. M. and Grest, G. S., Physical Review B fl, 4685 (1991-1).

96. Spencer, M. S., Review Article 525, 685 (1986).

97. Srolovitz, D. J. and Safran, S. A., J. Appl. Phys. 59, 247 (1986).

98. Srolovitz, D. J. and Safran, S. A., J. Appl. Phys. Q, 255 (1986).

99. Stauffer, D. Percolation Theory (1985).

100. Stein, J. and Krey, U., Z. PhysiltB 3_4, 287 (1979).

101. Stein, J. and Krey, U., Z. PhysiltB 51, 13 (1980).

102. Stillinger, F. H. and Weber, T.A., Phys. Rev. A 55, 978(1982)

103. Tam, S. W. and Johnson, C. E., J. Phys. A 59, L471 (1987).

104. Thompson, P. A., Grest, G. S. and Robbins, M. 0., Physical Review Letters 6_8, 3448
(1992)

105. Thouless, D. J ., Physical Review Letters 52, 1167 (1977).

106. Thouless, D. J. and Kirkpatrick, S., J. Phys. C 551, 235 (1981).

183

107. Umbach, C. C., Keeffe, M. E. and Blakely, J. M., J. Vac. Sci. Technol. A 2, 1014
(1991)

108. Weeks, J. D., in Ordering in Strongly F luctuating Condensed Matter Systems, edited
by T. Riste (Plenum, New York, N. Y.) 1980.

109. Weeks, J. 0, Physical Review B 25, 3998 (1982).

110. Wulff, G., Z. Kristallogr. Mineral, 5A , 449 (1901).

111. Yoshino, S. and Okazaki, M., Journal of The Physical Society of Japan 55, 415
(1977)

112. Yu, X., Duxbury, P. M., Jeffers, G. and Dubson, M. A., Phys. Rev. B 95, 13163
(1991-1).

113. Zangwill, A., Physics at Surfaces, Cambridge Univ. Press (1988)