.i .11 ..._ u. 1;.
.._.._..._3 3......“ . :

 

 

(.2 . , :5

z. . .{Léhul
w. .s .

..
52.13....
. .7)...‘
. 3...:

use}:
2. .V:I. .
. 3...... .
5:13:31)

.1. 3... I:
. rear... .5. x

2...... .
1.1, 1.

I...
:3...

1 x. 1... .3...
CF .. 11.5....
.21..)

. .3
.322:
. r. 1... x2...
33;...41)

21......
. :1 x

u ... ., {.3 i... t...
2:: Z I 3 .3
.
5;.

. .1}.
infra}...
.. 5.5....

.y. 1.1:.-. a

.11. :1: 2. .
4 . 3......

3.

u...£:(.......!..t
.7. a; to . . an
. a . . .2 ..

l
.5 .. .3
I...

5:
1...}
3:12,

P... 2...... z . y . 3.5.92
1:1:59w... 5 V x . . in...)
....r.

) 512...... :5
a» . a .2. £3.31. ,v
L3} .... .33....3. .. ,x.
.{..f..,.,..\ N)...1;.<.~
03‘ n)... .. xv

 

...(L .14 3...)...

«15...»... u

 

3:; n...
.f .3~ 1:).

 

 

.«v. 2....

 

I>~».v...h

L1“...

.51.... ..

, .....

 

i: E. .1.

. . . .tv‘txflnl.6u) .
9.... .1 x. i...
. . Hafiz... . i y

, L L

 

lllllNHIIllHlllllllllllllllllllllll{NHWINIIIINIINIIII

293 010206

This is to certify that the

dissertation entitled

STABILITY AND DYNAMICS
OF
LOW-fDIHENSIONAL SYSTEMS

presented by

SEONG GON KIM

has been accepted towards fulfillment
of the requirements for

PhD degree in Physics

mmy

P . M. Duxbury

Major professor

 

Date November 16, 1994

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

 

 

 

 

 

 

 

LIBRARY '
Michigan State
University

 

 

 

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

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

J

m
—]L 1|

MSU Is An Affirmative Action/Equal Opportunity inditution
Warns—9.1

 

 

 

 

 

 

i

 

 

 

 

 

STABILITY AND DYNAMICS
OF
LOW-DIMEN SION AL SYSTEMS

By

Seong Gon Kim

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

1994

 

ABSTRACT

STABILITY AND DYNAMICS
OF
LOW-DIMENSIONAL SYSTEMS

By

Seong Gon Kim

Understanding the structure and dynamics of physical systems has been one
of the major tasks of condensed matter physics. In this thesis, I study this issue
on four different physical systems: A superconductor connected to an external
current source, random networks under tensile stress, carbon fullerenes in contact

with an external heat bath and ferrofluid particles subject to an external field or

heat source.

The critical current density of the superconductor depends on the distribution
of its current density. For superconductors with defects, the critical current den-
sity will depend on how these defects enhance the local current densities. The
numerical simulation results as well as analytic predictions on the dependence of

critical current density on the geometry of the defects is presented.

The mechanical response of materials can be modeled by networks of nodes

and bonds. Tuning the energy functional of the bond deformation, I analyze the

 

strength of the system and simulate the fracture process. A scaling theory is also

presented.

Carbon is one of the most favored elements of scientists because of the many
remarkable properties it possesses. Discovery of carbon fullerenes spurred even
more interest in carbon clusters. I studied the thermal disintegration of fullerenes
within the framework of tight-binding molecular dynamics. Not only the melt-
ing temperature but also many intermediate structures including a previously

unknown “pretzel” phase have been identified.

A ferrofluid is a magnetic colloid made of magnetite particles suspended in
a liquid, usually a petroleum oil. When a magnetic field is applied, they form
well defined structures. Our study indicates that these particles show interesting
structures even in zero-field. We also find that the ring structure is very stable

and it should be observed in carefully controlled experiments.

 

 

 

 

 

To my dear wife
Jinwan
and my precious children

Danbee, Mara and Barum'e

 

 

 

 

ACKNOWLEDGEMENTS

I thank God for all the blessings He has stored for me and my family.

I would like to express my deep gratitude to my thesis advisors Professor
Phillip M. Duxbury and Professor David Tomének for their constant support, their
physical insights and for sharing with me their enthusiasm for physics. I especially
thank Dr. Duxbury for his generosity and pure concerns for my successful career
to allow me to work with Dr. Tomanek during his sabbatical leave. It was my
great honor and privilege to be an apprentice to such great teachers. During
the years of my learning to become a physicist myself, they have guided me with
confidence and love through the ups and downs of research. Their continual search
for simple and fundamental answers even in complex problems will always be my
guiding light.

I am also grateful to Professors S.D. Mahanti, W. Pratt, W. Repko, D. Stump
and W. Bauer for serving on my guidance committee. I would like to thank Mrs.
Janet King and Stephanie Holland for making my life in this department so much
simpler and enjoyable. The help from Mr. Ron Winsauer and Dr. George Perkins

in solving the computational problems also weighs greatly.

I wish to thank Professor J. Hetherington and Dr. Philippe Jund for their

numerous discussions, comments and help for the completion of the last chapter.

My deepest gratitude, however, can not go to anyone but to my dear wife
Jinwon and my precious three children, Danbee, Mara and Barunie. They have

been my strength, hope and the reason for completing this work. I sincerely desire

iv

 

 

 

that someday I can pay back even a tiny fraction of the countless sacrifices they

have made for me.

I also want to thank to all of my friends in this departments who made the
life in East Lansing so much more enjoyable and pleasant. My final thank goes to
Profs. W. Repko and W. Hartmann for keeping me in shape and recharging my

strength with weekend tennis workouts.

Financial support from the Natural Science Foundation and the Center for
Fundamental Materials Research of MSU without which this thesis could not

have been, is also greatly acknowledged.

 

 

Contents

Abstract ..................................

Acknowledgments ............................
LIST OF TABLES
LIST OF FIGURES

1 Introduction

1.1 Cracks and critical current ......................
1.2 Elasticity and fracture of disordered networks ...........
1.3 Melting of carbon fullerenes .....................
1.4 Stability of Magnetic Dipole Structures in Ferrofluids .......

2 Cracks and critical current

2.1 Introduction ..............................
2.2 Supercurrent enhancement at a crack tip ..............
2.2.1 The Model ...........................
2.2.2 Exact solution as a/A —» 0 ..................

iv

ix

xii

11

12

13

13

18

 

2.2.3 Saturation Effect as a/A —» oo ................

2.3 Numerical Simulation .........................
2.4 The Effects of a Crack on Critical Current .............
2.4.1 The critical state model ...................
2.4.2 Hotspot theory ........................
2.5 Conclusions ..............................

Elasticity and fracture of disordered networks

3.1 Introduction ..............................
3.2 Central force networks ........................
3.2.1 Effective medium theory for elasticity ............
3.2.2 Numerical simulation of damage and fracture .......
3.3 Scaling theory of tensile fracture stress ...............
3.4 Bond Bending Networks .......................
3.4.1 Background ..........................
3.4.2 Simulations of honeycomb networks .............
3.5 Conclusions ..............................

Melting of carbon fullerenes

4.1 Introduction ..............................
4.2 The tight-binding potential for carbon ...............
4.3 The recursion method ........................

22

25

32

34

40

40

42

43

45

48

57

57

60

69

71

71

76

84

 

4.4 Molecular dynamics .......................... 89

4.5 Simulation and Results ........................ 94
4.6 Conclusions .............................. 110
Stability of Magnetic Dipole Structures in Ferrofluids 111
5.1 Introduction .............................. 111
5.2 Modeling of the ferrofluid ...................... 114
5.3 Molecular dynamics of spherical tops ................ 122
5.4 Quaternion molecular dynamics ................... 126
5.5 Rings and chains ........................... 133
5.6 Molecular dynamics study ...................... 143
5.7 Ongoing work ............................. 147
5.8 Conclusions .............................. 152

 

List of Tables

4.1

4.2

5.1

5.2

Parameters for the band structure energy .............. 80
Parameters for the repulsive energy ................. 83
Units and conversion factors of physical quantities related to fer-

rofluids. ............... I ................. 119

Physical parameters used in ferrofluid study. ............ 120

 

 

List of Figures

2.1

2.2

2.3

2.4

2.5

2.6

2.7

2.8

2.9

3.1

3.2

3.3

3.4

3.5

3.6

A superconductor with elliptic defect ................. 15
Elliptic coordinates ........................... 17
The magnitude of current density near the crack tip. ....... 20
A superconducting lattice with line defect. ............. 26
The profile of the magnetic field near a crack. ........... 28
The profile of the magnitude of the current density near a crack. . 29

The current at the crack tip as function of crack size. ....... 30
The saturation current at the crack tip ................ 31
The Bean model (schematic). .................... 33
The stress-strain behavior of a single spring. ............ 46
The stress-strain behavior of a 50x50 triangular network. ..... 47
The fracture evolution (early stage) ................. 49
Intermediate configuration ....................... 50
The final fracture configuration. ................... 51
The Young’s modulus and tensile fracture stress .......... 54

 

3.7 Test of scaling of 0‘5( f) ......................... 55
3.8 Cu and p. for a honeycomb tube. .................. 61
3.9 The stress-strain curve for honeycomb tube. ............ 62
3.10 The initial configuration of tensile fracture of a hexagonal lattice . 64
3.11 The final path of tensile fracture of a hexagonal lattice ...... 65
3.12 The response of a graphite sheet containing a crack-like defect cluster. 66

3.13 The average fracture stress as a function of void volume fraction . 67

4.1 Carbon fullerenes ............................ 72

4.2 Comparion of recursion method vs. conventional tight-binding

method. ................................ 90

4.3 Temperature dependence of the total energy E per atom for the
C60 fullerene. ............................. 96

4.4 Specific heat cv(T) per atom of 020, C60 and 0240 .......... 97

4.5 “Snap shots” illustrating the geometry of a 060 cluster at temper-

atures corresponding to the different “phases” ............ 99

4.6 Atomic binding energy distribution in a 060 cluster at temperatures

characteristic of the different “phases”. ............... 100

4.7 Bond length distribution in a 060 cluster at temperatures charac-

teristic of the different “phases” ..... ' ............... 102

4.8 Bond angle distribution in a Cso cluster at temperatures character-

istic of the different “phases”. .................... 103

 

 

 

4.9 Mass spectrum in a C50 molecule at temperatures characteristic of

the different “phases”. ........................ 104

4.10 The entropy 3 per atom for the 060 fullerene and the metastable

Cso chain. ............................... 108
5.1 Soft-core potential for a ferrofluid particle. ............. 117
5.2 Two dipoles on :r-axis. ........................ 121
5.3 Definition of Euler angles. ...................... 123
5.4 A chain and ring of 10 magnetic dipoles. .............. 134

5.5 Energy of a chain and a ring of 10 magnetic dipoles in the absence

of external field ............................. 135
5.6 A vector representation of a ring of 10 magnetic dipoles. ..... 137
5.7 Rings of 20 magnetic dipoles in the external field. ......... 140

5.8 Energy of a chain and a ring of 10 magnetic dipoles in the external

field. .................................. 141
5.9 Lower critical field for magnetic dipole rings ............. 144
5.10 Upper critical field for magnetic dipole rings ............. 145
5.11 Energy barrier E' of ferrofluid rings. ................ 148
5.12 Estimation of melting temperature T5, of ferrofluid rings ...... 149
5.13 Aggregation of 64 magnetic dipoles .................. 151

 

 

 

 

Chapter 1

Introduction

The major task of theoretical physics is to understand the fundamental laws which
govern nature, use them to explain observed phenomena and predict the properties
of new physical systems. The subjects of interest for theoretical condensed matter
physics are electronic, magnetic, structural and dynamical properties of solids and
liquids. The most fundamental technique to address most of these problems from
first principles is the ab initio Density Functional Theory (DFT) [Lund83] in the
Local Density Approximation (LDA) [Hohe64, Kohn65], which has been shown to

describe the ground state properties of solids with high accuracy.

The heavy computational requirement associated with ab initio methods such
as the LDA, however, limits ‘the range of its applicability to relatively small sys-
tems with a high symmetry. In order to describe the dynamical behavior of large
systems with low symmetry or disordered systems, we need to use somewhat less
demanding approximation schemes. In this thesis, I use the Tight-Binding (TB)
method, or parameterized Linear Combination of Atomic Orbitals (LCAO) for-
malism to study the dynamics and stability of carbon fullerenes. TB method

conveniently parametrizes ab initio LDA results for structures as different as Cg,

 

 

 

carbon chains, graphite, and diamond [Toma91].

This method is much more efficient in terms of the computational effort and has
the essential quantum mechanical features still in the formalism. TB method can
be made yet more efficient and fast by adopting the recursion technique [Zhon93]
in diagonalizing the Hamiltonian matrix. With an exact numerical diagonalization
of the Hamiltonian matrix, the CPU requirement scales as N 3, for N the number
of particles. Using the recursive method makes the computational load scale as
N. This factor makes a huge difference when we put more than a few hundred
atoms in the system. Furthermore if you need to do a Molecular Dynamics (MD)
study of these systems the benefit of a linearly scaling numerical method will let

us afford to simulate large systems for a sufficiently long time span.

When the object of our study becomes bulk materials, such as a piece of
superconductor connected to an external current source or a piece of metal bar
clamped in a tensile-test machine, however, there are too many particles in the
system to be analyzed or simulated by even the TB method. Moreover, when
these bulk materials have random defects it would be hopeless to treat them as the
collection of individual atoms or molecules. For systems like a dielectric medium
or a superconductor, I use the continuous medium approximation and simply solve
the relevant equations of state parameters such as Maxwell’s equations and/ or the

London equation.

For mechanical stability calculations, many systems can be modeled by net-
works of nodes and bonds. A node denotes a point fixed in the bulk material and
a bond represents the interaction between these nodes. By postulating the energy

functional for these bonds and tuning the parameters therein this network model

 

 

 

can reproduce the mechanical response of the bulk materials reasonably well. In

the studies included in this thesis, I use one of the simplest energy functions, the

“simple harmonic oscillator” model. The bond between nodes is represented by a

linear Hooke’s spring. The strength of the material is represented by the strength
of the springs and the defect can be introduced easily by modifying individual
spring constants. The external stress or deformation is specified by setting the
prOper boundary conditions. The optimum configuration of the system is obtained
by minimizing the total energy of the network using the Conjugate-Gradient (CG)
method. For this method to be successful, the number of nodes should be suf-
ficiently big and it is usually a couple of orders of magnitude bigger than the

number of atoms we can deal with in tight-binding MD or similar methods.

During my years of training for completion of the Doctoral degree, I have
worked on four major different subjects. Each subject is in a separate chapter
and each chapter has its own introduction and conclusion. In the following, I

introduce each chapter and its subject briefly.

1.1 Cracks and critical current

In Chapter 2, the effect of defects on the critical current density of superconductors

is studied.

After Kamerlingh Onnes reported the first observation of superconductivity in
April of 1911, superconductivity became one of the most researched subjects in
the science community. There has been vast developments in terms of discover-
ing more and more superconducting materials and enhancing the quality of the

superconductors. It is a well known fact that superconducting materials are not

 

 

 

 

 

 

always superconducting. They need to have very special conditions to be super-
conducting. One of these conditions is the so-called “critical current density”,
jc, of a superconductor. There is a maximum current density a superconductor
can support. Thus when you have an inhomogeneous current density distribution
in the superconducting medium, the location where it has the maximum current
density is the most vulnerable to this transition. It is well known that cracks (of

width, w > £,) are effective in enhancing local current density.

In this chapter I calculate the size of these hotspot supercurrents as a function
of the crack length, a, and the London penetration depth, A, using 2D London
theory. The current density enhancing factor for a given crack size is studied as a
function of the penetration depth and the intensity of vortex pinning. We argue
that large local supercurrents near a surface crack nucleate vortex creation. If
flux pinning is weak enough these vortices flow under the influence of the large
local Lorentz force. The dissipation so produced can lead to a reduction in ob-
served critical current. If flux pinning is moderate, the first additional vortices
nucleated near the crack tip are pinned, in a region we label p, the flux pinning
zone. In the case of a through crack in a thin film, we then argue that jc(a)
reduces as jc(a)/jc(0) ~ (p/a)‘ (for a/L << 1), where L is the film width and
z = 1/2 for the simplest London theory. We compare this theory with the Bean
(critical state) model which predicts that the critical current is (approximately)
determined by the cross-section available for supercurrent, so that in a film con-
taining a crack, jc(a)/jc(0) ~ 1 — a/ L (ignoring self field effects). We argue that
superconductors with sufficiently weak pinning should obey the hotspot theory,
while sufficiently hard superconductors should obey the critical state model, and

suggest experiments that should illustrate these two limiting cases.

4

 

 

 

 

 

1.2 Elasticity and fracture of disordered net-
works

In Chapter 3 I discuss the elasticity and mechanical failure of random spring

networks.

Many alloys and composites are disordered. In fact no material is perfect and
the strength of a given material is very much dependent on these minor flaws
or defects. Thus understanding the effects of random defects and their statistics
is very crucial in analyzing material strength. The effective elastic properties
of these materials can be calculated using effective medium theory and other
“homogenization” methods. Due to the fact that fracture is initiated in regions of
high stress, it is to be expected that homogenization methods, unless they focus
in the crack growth region, should be poor predictors of the effect of disorder on
failure strength.

Random networks with “realistic” interatomic and many body potentials
are now routinely simulated on a computer. Recently, detailed analytic and
numerical analysis of idealised disordered microstructures have been performed
[Herr90, Sahi9l]. Along with new scaling theories specifically taking into account
the extreme statistics of brittle fracture[Duxb87a, Duxb91, Bea188], these simula-
tions have lead to new insights into the fracture of materials. In this chapter, we

discuss the use of these methods to study fracture of brittle materials.

 

T—_——i

1.3 Melting of carbon fullerenes

Chapter 4 contains the work I did for carbon fullerenes. Carbon fullerenes, such
as Can [Krot85], are intriguing molecules with a hollow soccer ball structure and a
large range of interesting properties [Bill93]. Due to the high structural stability of
these systems [Cur188], almost no information is available about their equilibrium
phases at temperatures close to and exceeding the melting point of graphite.
We expect the phase diagram of these strongly correlated structures to be very

interesting and to consist of several well identifiable “molten” phases.

A second reason to study the high-temperature behavior of fullerenes is re.
lated to their formation. While bulk quantities of fullerenes can now be rou-
tinely synthesized in a carbon arc [Krat90], the microscopic mechanism of their
aggregation from the gas phase is still the subject of a significant controversy
[Waka92, Waka93, Ebbe92, Kern92, Held93, Hunt94]. During the formation pro-
cess, fullerenes will be created from many different intermediate structures. They
also will have many unsuccessful trials and to form such a high symmetry struc-
ture from the gas phase will require a tremendous amount of encounters between
even correctly matching parts of the structures. Therefore to simulate the forma-
tion of fullerenes we will have to have many carbon atoms in a simulation box
and wait for a long time to have some closed structures, without even counting
the number of atoms in the cluster. As an alternative approach, we hope that a
detailed study of the annealing process may shed new light on this controversy,

since intermediate structures, which occur during thermal quenching in the rare

 

gas atmosphere, may also be observed while heating up these structures to high

temperatures.

 

 

 

 

A third motivating reason for this study are recent collision experiments which
indicate that fullerene molecules are extremely resilient and only fragment at en-
ergies exceeding as 200 eV [Camp93]. Since in the collision process the kinetic
energy is transferred into internal degrees of freedom as “heat”, additional infor—
mation about the intermediate structures of colliding fullerenes could be obtained

by investigating those of the superheated molecules.

In an attempt to elucidate the above three problem areas, we performed a
detailed molecular dynamics (MD) study of the melting and evaporation pro-
cess of three prototype fullerenes, namely 020, 060, and 0240. In our simula-
tion, we investigated the response of a canonical ensemble of fullerenes to grad-
ually increasing heat bath temperatures using Nose-Hoover molecular dynamics
[Nose84, Hoov85, Alle90], complementing recent microcanonical ensemble simu-
lations on 060 and Cm [KimE93]. Of course, the quality of simulation results
depends primarily on the adequacy of the total energy formalism applied to the
fullerenes. Since our simulation also addresses the fragmentation under extreme
conditions, the formalism used must accurately describe the relative stability of
fullerenes with a graphitic structure and their fragments, mostly short carbon
chains. The corresponding transition from .31:2 to sp bonding can only be reliably
described by a Hamiltonian which explicitly addresses the hybridization between
02.9 and 02p orbitals. A precise evaluation of interatomic forces is required to
obtain the correct long-time behavior of the system. On the other hand, the long
simulation runs necessary to equilibrate the small systems during the annealing
process, and the large ensemble averages needed for a good statistics, virtually

exclude ab initio techniques as viable candidates for the calculation of atomic

forces.

 

 

 

 

For this reason, we base our force calculation on a Linear Combination of.
Atomic Orbitals (LCAO) formalism which conveniently parametrizes ab initio
Local Density Approximation (LDA) [Hohe64, Kohn65] results for structures as
different as Cg, carbon chains, graphite, and diamond [Toma91]. We use an adap-
tation of this formalism to very large systems, which is based on the fact that forces
depend most strongly on the local atomic environment and which has been imple-
mented using the recursion technique [Zhon93]. The force calculation can be per-
formed analytically to a large degree. This not only speeds up the computations
significantly, but also provides an excellent energy conservation AE / E 5.. 10‘10
between time steps in microcanonical ensembles, and a linear scaling of the com-
puter time with the system size. Most importantly, forces on distant atoms can

be efficiently calculated on separate processors of a massively parallel computer.

Using MD simulation, we heat up carbon fullerenes and monitor the changes
in various structural and thermodynamic responses of the system. We identified
many well defined intermediate structures of fullerenes. Among others, we found
a most dramatic transition to a previously unknown “pretzel” phase to occur at
a high temperature of 132.4000 K. The temperature where fullerenes disintegrate

into carbon chains is explained using quantitative results for the entr0py.

1.4 Stability of Magnetic Dipole Structures in
Ferrofluids

Chapter 5 contains the work on ferrofluids. A ferrofluid is a magnetic colloid
made of magnetite particles suspended in a liquid, usually a petroleum oil [Rose85,

Wang94]. A typical electrorheological (ER) fluid consists of a suspension of fine

 

 

 

 

dielectric particles in a liquid of low dielectric constant [TaoR92, Bloc87, Fili88].
The aggregation of these particles has been the subjects of many experiments
and numerical simulations. Experiments find that upon application of an electric
field, dielectric particles in ER fluids rapidly form chains which then aggregate
into thick columns [Hals90, Chen92, Mart92, Gind92]. Magnetic ferrofluids forms
a Similar structure of chains and thick columns. The dependence of the periodicity
of columns on the sample thickness has been recently studied using the magnetite
colloids [Wang94]. Also many numerical simulations using the molecular dynamics
(MD) [Ha5389, Klin89, Kusa90, Bonn92, TaoR94] has been carried out to find the

orientational ordering in ER liquid crystals.

Even though a few MD simulation results have been published for ferroelectric
dipole liquids [WeiD92a, WeiD92b], almost no MD simulation results are available
for the magnetic ferrofluid systems except a few simulations using the Monte Carlo
technique[Weis93, Stev94]. There is one fundamental difference between ER fluids
and magnetic ferrofluids even though they have many common properties. In ER
fluids, the dipoles are induced by an applied electric field and they are always along
the field. The strength of their dipole moments depends on the local electric field
and therefore they do not show any interesting properties in the absence of the
external field. On the contrary, the magnetic particles in ferrofluids have their
own pre-assigned magnetic moments. This means that in general they are not
obliged to be aligned with the external field and we have to consider the rotations
of these particles as well as their positions. Also we should expect that these fluids

have some non-trivial structures even in the absence of the external field.

In this chapter I will discuss the early stages of colloidal aggregation in fer-

rofluids and the stability of the ring structure of these particles. We found that
9

—¥—

 

 

 

the ring is more stable for ferrofluid particles in low field and low temperature as
long as it has more than 4 particles per ring. An estimation of the energy barrier
between the ring structure and chain structure is presented along with the melting

temperature in an external magnetic field.

10

 

 

 

 

Chapter 2

Cracks and critical current

In superconductors, cracks (of width, 1.0 > 5,) are effective in enhancing local
supercurrent, and we calculate the size of these supercurrent hotspots as a function
of the crack length, a, and the London penetration depth, A, using 2D London
theory. In the A —-> oo limit and constant injected current density, we show
that this 2D solution is also ezact for the current flow in a thin film containing
a through crack. We argue that large local supercurrents near a surface crack
nucleate vortex creation. If flux pinning is weak enough these vortices flow under
the influence of the large local Lorentz force. The dissipation so produced can

lead to a reduction in observed critical current. If flux pinning is moderate, the

 

first additional vortices nucleated near the crack tip are pinned, in a region we

 

label p, the flux pinning zone. In the case of a through crack in a thin film, we
then argue that jc(a) reduces as jc(a)/jc(0) ~ (p/a)‘ (for a/L << 1), where L
is the film width and a: = 1/ 2 for the simplest London theory. We compare this
theory with the Bean (critical state) model which predicts that the critical current
is (approximately) determined by the cross-section available for supercurrent, so
that in a film containing a crack, jc(a)/jc(0) ~ 1 —a/ L (ignoring self field effects).
We argue that superconductors with sufliciently weak pinning should obey the
11

 

 

 

hotspot theory, while sufliciently hard superconductors should obey the critical state

model, and suggest experiments that should illustrate these two limiting cases.

This chapter contains materials which have appeared in the following publica-

tion:

a [KimSQl] S.G. Kim and RM. Duxbury, Cracks and critical current, J. Appl.
Phys. 70,3164 (1991).

2.1 Introduction

In the Bean[Bean64] (critical state) model of hard type II superconductors, the
supercurrent flows in shells in which the current is at its critical current level, and
the remaining parts of the superconductor carry no current. The macroscopic
critical current occurs when the “local” critical current flows throughout a cross—
section of the sample. This model predicts that the observed critical current is
relatively insensitive to extended defects such as cracks (critical current is reduced
roughly according to the cross-section of the extended defect). In this paper we
show that if pinning is sufficiently weak, a crack can cause a much more rapid
reduction in observed critical current due to the nucleation and flow of vortices

from the crack tip.

We initiate a study of the origins and consequences of large local supercur-
rents, by studying the effect of a crack on supercurrent flow in thin film super-
conductors. This is a natural starting point, as it is known from the study of
electrical[Duxb87a, Duxb87b] and mechanical failure[Kell86] that cracks or crack
shaped flaws are very efficient at enhancing current or stress. As the simplest non-

trivial illustration of the effect, we use 2D London theory to show that a crack

12

 

 

 

can lead to large local supercurrents. This effect is demonstrated analytically in
Sec. [2.2] and the analytic calculations are supported by numerical solutions to
the London equation in Sec. [2.3]. In the limit A H 00, we show that the 2D

London theory is exact for a thin film containing a through crack.

In Sec. [2.4] we consider the effect of large local supercurrents on critical cur-
rent. We argue that large local supercurrents act as nucleation sites for vortex
creation. Once created, the vortex lines will flow and cause dissipation, if the local
flux pinning is insufficient to resist the Lorentz force on the vortex line. We thus
argue that large local supercurrents degrade critical current most in weak pinning
superconductors where flux lines can be depinned by the additional Lorentz force
due to large local supercurrents. Using these ideas in combination with the Lon-
don theory calculations of Sec. [2.2], we make some specific predictions for the
dependence of the critical current of a superconducting film on the length of a
(lithographically drawn) through crack in the film. These predictions are com-
pared with the results one would expect using the simplest critical state model.

Sec. [2.5] contains a summary of the main results of this study.

2.2 Supercurrent enhancement at a crack tip

2.2.1 The Model

To illustrate the effect of crack on supercurrent flow, we use 2D London the-
ory. London theory is known to provide a good semi-quantitative description of
strongly type II superconductors away from Te [Tink65, Duze81, Ti1186, Orla91].
In particular the London model does a fairly good job of predicting the field profile

and current flow near vortex lines. A crack can be modelled as a line of vortices,

13

 

 

 

 

 

 

 

 

 

so we expect that London theory will also provide a good first approximation in
this problem. We must keep in mind though that the theory is incorrect on length

scales shorter than 5..

Although the 2D theory is only strictly valid for an infinitely thick film con-
taining a through crack, we show later in this section, that in the A —+ co limit, it
is also exact for the current flow pattern in a thin film containing a through crack.
In the 2D geometry, current flow is in the z-y plane, so that only the z-component
of magnetic field is finite. This reduces the London theory to a scalar Helmholz

equation for the z-component of the magnetic field, B:
V2805?) = 3(1” yl/Az' (2'1)
Here B is a function of :e and y, and the current is found from B via the Maxwell

equation,

C

4WV x [0,0,B(z,y)]. (2.2)

j:

Using Eqs. (2.1) and (2.2), we will calculate the current flow in the geometry

shown in Fig. 2.1.
We wish to calculate the current flow pattern using as a boundary condition
lbw) = loijoexp(—z/’\)l as y "i too' (2‘3)
This is achieved by imposing the field boundary condition,
_ B(z,y) = 30(3) = Aoexp(-:r/A) as y —» ioo. (2.4)
with
41Ajo

A0 = . (2.5)
14

 

 

 

 

 

 

 

Figure 2.1: A superconductor with elliptic defect used in the analytic calculations.

15

 

 

 

 

 

 

 

Another boundary condition we impose is that there is no current normal to the
surfaces of the superconductor (including the surfaces of the crack-like ellipse).
Note that since we are solving for a problem with constant current boundary
conditions inside the film and zero current outside the film, we in principal solve
for the current pattern inside the film first and then generate the magnetic field
patterns interior and exterior from it. It turns out that it is mathematically
convenient to still work in terms of the magnetic field, but the physical geometry
should be kept in mind when thinking about the solution. To satisfy the boundary
conditions at the crack surface, it is convenient to use to elliptic coordinates

(£2flllMors53] as illustrated in Fig. 2.2. In this coordinates an ellipse defined by:
2:2/a2 + yz/b2 = 1 (2.6)

using the usual elliptic variables,

a: = Ccoshécosn (2.7)
y = C sinhfsinn

where
a = C cosh 60
b = C sinh £0 (2.8)
0 = m

Physically, f and n are the radial and angular coordinates respectively and the
surface of an ellipse is generated by fixing f = £0 and allowing 1) to vary between

—1r/2 and 1r/2.

Separation of variables:

B(£,n) = f(£)g(n) (2.9)
16

 

 

Figure 2.2: Elliptic coordinates.

leads to

22.9.

dnz + (s + 2q cos 21])g = 0 (2.10)
and

at2

32% — (s + 2q cosh 2£)f = 0 (2°11)

where q = (C'/2A)2 and s is a separation constant. Equations (2.10) and (2.11)
are the original and modified Mathieu equations[McLa46, Abro72, Grad80] and
arise in a variety of quantum and classical scattering problems[Bowm87, McLa47]
(note though that the sign of q is opposite that occurring in such problems). To

include the boundary condition of Eq. (2.4), we write

B(£,n) = 30(3) + 3105,17) (2.12)

and we choose radial solutions to Eq. (2.11) that ensure that 31 (6,1)) goes to zero

at large 5.
2.2.2 Exact solution as a/A —-» 0

Equations (2.10) and (2.11) do not yield to simple closed form solution, so it is
instructive to consider the limit a / A -> 0 where a closed form solution is straight-

forward (this limit is also relevant to thin films where the effective penetration

depth is large). In this limit, the boundary condition
30(2) = Ao(l - z/A) for y -i :too (2.13)
leads to

B(z,y) = 1400-3”)

+Ao(a//\) exp(-[£ - fol) cos 17 (2-14)
18

 

:2"

 

 

 

 

for the magnetic field exterior to the ellipse. From this solution it is straightfor-

ward to calculate the supercurrent pattern in the film using Eq. (2.2). The results

je(£,n) = (fa/THC coshtsin 1) - aexpl-(é — (0)] sin 17} (2-15)
and

ME, 17) = (fie/THC sinh t cos n + aexpl-(t - £01608 71} (2-16)
where

r2 = 02(sinh2£ + sinzn). (2.17)

A detailed study of these solutions near the crack tip shows that there is a large
enhancement in the current density in the y direction at the crack tip. The current
density may be calculated as a function of distance (in the z-direction - see Fig.
2.1) from the crack tip (on the central axis of the crack) by evaluating the current
in the 1) direction as a function of f: Using f = cosh-1(z/C') at n = 0 and taking
a: = a+r, we find the current behavior presented in Fig. 2.3 (the result is identical

to that found in a classical conductor containing an insulating inclusion[LiYS89]).

The behavior of Fig. 2.3 can be asymptotically summarized as follows:

1 + (a/2p)1/2 if r < p;
ju(")/j0 ~ 1+ (er/2r)”2 if p < r < a; (2-18)
1+a2/(2r2) ifr >> a;

 

 

 

4)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

10 TTTT TT—T—l— r1. Tr I l l l l I I—
10 L r r I I 1
F 1
105 — -—
.1
r .
o - .
”-5 ..
\ 10° L ——
A r- 1
i: . 1
>5 r- -(
H 10‘5 -- —
L * .1
)- a
- -4
r III
’4 1 14 L L 1 1 LLI 1 Li 1 1 1 1 l 41 4

10"10 10"5 10° 105 1010

1"

Figure 2.3: The magnitude of current density as a function of distance (in the
z-direction) from the crack tip when A > a and a/b = 1000. The asymptotic
behaviours predicted by the analytic calculation of Eqs. (2.18) are clearly evident.

20

 

 

 

where p = b2 / 2a is the crack tip curvature. One important feature of this result

is that the magnitude of current in a region of size p near the crack tip is

 

_I_t_ip_ : lop jy(r)dr
Io f0” J'o exp (-2=/ AW
~ 1 + (a/2p)1/'-’. (2.19)

The crack tip current increases as the square root of the crack length and becomes

very large for large a.

Now consider a film of finite thickness in the z-direction, which has the same
cross-section (that of Fig. 2.1) for every value of z. The current solution found
above solves this problem also, as there is no z-direction current, and the boundary
conditions in the z-y plane are the same as for the 2D case. Every plane of the
film thus has the same current pattern as the 2D problem, and there is no :-
dependence to the magnetic field pattern inside the film. Of course outside the
film there is a very strong z-dependence in the magnetic field generated by the
current sources in the film. The magnetic field exterior to the film thus has z,y,
and 2 components that depend on 2,3] and 2. We are interested in the effect of
the current pattern inside the film on critical current, and will concentrate on the

interior solution for the remainder of this chapter.

For finite A, the interior current pattern problem is truly three dimensional.
We also expect the current enhancement effect [see Eqn. (2.19)] to be reduced,
as current is already concentrated in a region of thickness A near the surface of
the superconductor. In fact, the reader may intuitively expect that for a >> A,
1.5,, should saturate at large a. This intuition turns out to be correct, and we
have quantified this “saturation effect” by developing asymptotic series for 1,3,, as

a/A -> 00. We use 2D London theory, and so that the results we will derive are

21

 

 

 

only rigorously valid for infinitely thick films. We do however expect that they
will also provide a good first approximation for thin films when A is large, as we
know from the above exact solution, that as A -—> co the 2D solution becomes

exact .

2.2.3 Saturation Effect as a/A —-> 00

In this limit, the angular function is rapidly decaying in 17. In addition, we are
interested in the limit b/a —> 0 so that 60 = tanh‘1(b/a) ~ b/a is small. It thus
makes sense to expand Eqs. (2.10) and (2.11) for small 1) and 5 and to use these
expansions to develop approximate solutions for the magnetic field at the crack

tip. A second order expansion in both 6 and 7] leads to

3217—3- + (k — mm = 0 (2-20)
and

g — (I: + 135’) f = 0. (2-21)
Choosing

I: = (2n + 1)]: and h = 0/1 (2.22)

we find that g and f are given by

9(a) = exr(-hn’/2)Hzn(\/7m) (2.23)
f(£) = Y(2n+1/2,~/2_h£) (2-24)

where H..(:c) is a Hermite polynomial of order n, and Y(n, z) is a parabolic cylinder
function[Abro72, Grad80] of order n [We have chosen the even Hermite polyno-

mials to ensure the correct symmetry about 1) = 0]. Both of these solutions can

22

 

 

 

 

 

be chosen to be rapidly decaying, which suggests that they will provide a useful
representation of the true solutions near the crack tip. (Further mathematical
justification for the truncation leading to Equations (2.20) and (2.21) may be
obtained from a Sturm-Liouville analysis[Mors53].) Using these functions as our

basis set and satisfying the boundary conditions, we find,

8(8) ll) ~ 30(2)

Moi3..exp(-hn2/2)H2.(¢in)r(4"; 1.¢2’h:) (2.25)

 

n=0

with

 

£22"(2n)!v(4"2+ 1,\/'27{eo)s,. =

/: exp(-hnz/2)Hzn(\/l;n)[l — exp(—§ cos 11)]. (2.26)

Here we have extended the angular integration to infinity with negligible error, as
Hn(n) is rapidly decaying in 1]. In the limit a/ A —-> co, the integral on the right
hand side of Eq. (2.26) reduces to

[0° eXP(-hn’/2)Hzn(\/hn)dn = 3,1591%! (2-27)

We then find that for a/A —2 00,

B(£.n) ~ 30(3)
a. 2 exp(_hq2/2)Hzn(\/in)Y[(4n + 1V2, flit] (2.28)

 

 

+A° ”2:; 22"n!Y[(4n + 1) / 2, V5360)
and hence
MM = 0) ~ 2'06““
, °° Y’[(4n + 1)/23 MC] (229)
+J° 2;, E"Y[(4n + 1)/2, J2—hto]

23

 

where

fi(—1)n(2n)1

with
Z 6,, = 1. (2.31)
n=0

Exactly at the crack tip, the series (2.30) is very slowly convergent (~ 11‘3/2
at large n) and 100,000 terms are required to achieve 3 figure accuracy for B
right at the crack tip. For any finite distance, p from the crack tip, however,
the series (2.28) is slowly convergent only up to terms N " ~ A/ p, after it is
exponentially convergent. In a similar way, the current density series (2.29) is
poorly convergent (marginally convergent at best) right at the crack tip, but is
exponentially convergent for large 12 provided p is finite. In applying our results
to experiment, or in comparing with the numerical work described in Sec. [2.3],
we must introduce a cut-off or coarse graining into the continuum theory. A study

of the series (2.29) integrated up to a cut-off or lattice spacing A yields
12 N Baa,» = 0) — 812,17 = o) (2.3,)

Io Bo(z = 0) —- Bo(z = A)

where {0 ~ b/a and {c = cosh’1 (a + A) / C . Since A > A, the parabolic cylinder

function may be approximated by its asymptotic behavior. We thus find,

Itip _ Isat
73w) _. co) = 10

1 — 23;. e.exp{—(/(4n + IMHO + 2a)“ — 1]}
'” 1 «xx—zaps) '

 

(2.33)

24

 

Here we have written A = Bp, and fl is the only free parameter in comparison
with the numerical simulations of the next section. Finally, the scaling behavior

of Eq. (2.33) for p < A < a is found to be

’7': ~ (AW/211 + moss/p) (234)

where 63 is a constant of order 1.

2.3 Numerical Simulation

We have used numerical simulations to test the analytic theory of the previous

section, and to calculate the current and field profiles around the crack tip. We use

a L X L square lattice, finite difference representation of the Laplacian operator,

and put the two line defects of length a lattice spacings and width 1 lattice spacing,
symmetrically at the both sides of the 2-dimensional array (see Fig. 2.4). The
symmetry of the system allows us to work with a system 1 / 2 the size of the actual
one. We inject current into the system from the bottom in the y direction. When
there is no defect the field and current profiles can be solved analytically and the
solutions are

sinhuL — 21/11

B = A° sinh(L/A)

, (2-35)

j = “might/3m 17- (2.36)

To simulate the effect of a crack, we ensure that the magnetic field a long
Way from the crack is that of Eqn. (2.35). We also fix the magnetic field at
surface lattice sites and at the surfaces of the crack. This ensures that no current

25

 

 

Figure 2.4: A superconducting square lattice with line defect used in the numerical

calculations.

26

 

l w.

 

flows across these boundaries. With these boundary conditions, the magnetic field
profile is calculated using the conjugate gradient method (For more details about
the method, see Ref. [Pre586]).

The square lattice for our numerical simulation is depicted in Fig. 2.4 and
the numerical result showing the magnetic field profile and magnitude of the
current density are plotted in figures 2.5 and 2.6. The magnetic field profile has
an exponential decay away from the crack tip (although it is faster than at a free
surface). The current, however shows a marked peak near the crack tip. This
current crowding near the crack tip is the effect of primary interest here, and
was calculated analytically in the previous section, Eqs. (2.18), (2.19), (2.33) and
(2.34) [Current is calculated from the difference of the magnetic field at adjacent
sites - the lattice version of Eq. (2.2)]. We plot the current at the crack tip for

various penetration depths and crack sizes in Fig. 2.7.

In this figure, the current is normalized to that occurring in a bond at the
free surface a long way from the crack. For large A, and a < L, the current
enhancement obeys the square root increase of Eq. (2.19). When a gets close in
size to L, finite lattice effects lead to a sharp rise in the current enhancement. For
smaller A, the current enhancement saturates for A < a < L, and again shows a
sharp increase as a approaches L due to finite size effects. The value at which the

current enhancement saturates (for a < L) is presented in Fig. 2.8.

In this figure, we also plot the predictions of Eq. (2.34) with ,6 = 1. It is seen
that the theory of the previous section provides an excellent fit to the numerical
data, and this justifies the use of the equations (2.19) and (2.34) in the scaling

analysis of critical current presented in the next section.

27

 

 

) 1.

 

 

 

 

\ 1
\\‘\‘....III I’
I,

\\ \\ \‘1..v’l
\ \ 1 a
(1 \ with,"

( \\ \‘ .10
\\\“

 

Figure 2.5: The profile of the magnetic field near a crack, found using numerical
simulations on 200 x 120 square lattice array for the case of A/ A = 10, a/ A = 30
and b/ A = 1.

28

 

 

Figure 2.6: The profile of the magnitude of the current density near a crack, found
in same numerical simulations for Fig. 2.5.

29

 

7-, 1111-11.. _-_.-- _.

 

 

 

 

 

 

 

 

10] T ED
9.
5..
7r-
6)-
5)-
4..
.93»

\D‘ O

.332- 03

o §°

0 g x x

. X

x a D D

If 1 L l
l

Defect Size (a/A)

Figure 2.7: The current at the crack tip (Lip/Io) as the function of the crack size
(a/A) for : A/A = oo(®),50(+),10(0),4(><),2(Cl).

30

 

 

 

 

 

100 1 7 WIFIIT I WWI—Tm

11'71

50

I
111111

I
I

 

 

 

 

1 5 10 so 100

Figure 2.8: The saturation current (1.») at the crack tip (see text) when A <
a < L. The solid line is obtained using Eq. (2.33).

31

 

 

 

 

 

 

 

2.4 The Effects of a Crack on Critical Current

In this section, we use the theory of Sections [2.2] and [2.3] to predict the effect
of a crack on the observed critical current of a superconducting film. In all cases,
we consider current injected into the film in a direction perpendicular to the long
axis of a crack (see Fig. 2.9), that extends all the way through the film (a through
crack). We first consider the prediction that the simplest critical state model
makes for this geometry. We then go on to a discussion of the critical current
reduction to be expected, if the large local supercurrent at the crack tip initiates

the dominant dissipative mechanism.

2.4.1 The critical state model

In the critical state model[Bean62, Bean64] of hard superconductors, the critical
current occurs when the current density is at its critical value all the way across
the sample. At lower values of external field or current, the critical current density
only flows in a shell near the surfaces of the film as depicted for a 2D geometry
in Fig. 2.9.

As the external current is increased, the current shells increase in size until they
fill the smaller cross-section extending outwards from the crack tip. Still ignoring
self-field corrections, the film case is also straightforward. Current shells now occur
on all edges of the film, and there is a z-dependence to the magnetic field and
current flow, but the total current carrying capacity of the film is still proportional

to the smallest crossection available for superflow. The critical current reduction

 

 

 

 

 

 

 

 

 

 

 

 

 

 

r INPUT CURRENT

 

Figure 2.9: The Bean model (schematic). Current flows in the shaded region.
Here we have rounded the sharp edges inherent in the model on a length scale of
A, as is expected on physical grounds.

33

 

 

 

due to a through crack in the geometry of Fig. 2.9 is then given by

%l ~ 9—2—9 (2.37)

ignoring self-field effects. This result shows a very weak dependence on crack
length compared to that of the hotspot theory discussed below. Including self-field
corrections (for example within the Kim model[KimY62, Ande62]) does modify
the result (2.37) and makes the supercurrent patterns more complex [Bean64] but

the crack sensitivity of the model is still weak.

2.4.2 Hotspot theory

We wish to put forward an alternative perspective on the effect of a crack on
critical current. Consider the limit of weak pinning, and slowly increase the current
injected into the film containing the crack (again the direction of current flow is
perpendicular to the long axis of the crack). Due to the current enhancement
effect at the crack tip [see Eqns. (2.18) and (2.34)] vortex nucleation will first
occur at the crack tip. If flux pinning is weak, the vortices so nucleated will flow
causing dissipation. (This has been explicitly seen in simulations of current flow
in the closely related problem of a crack in a 2D Josephson junction containing
a crack [XiaW89]). As pinning increases, the first vortices nucleated at the crack
tip will be pinned, and the film will support a higher external current before
dissipation sets in. As pinning increases further, a larger number of the newly
nucleated vortices remain pinned in the nucleation zone near the crack tip. We
label the zone in which these vortices reside the pinning zone, and ascribe to it
a radius p (note that we assume that there is no difference between the pinning
strength of the superconductor at the crack tip compared to elsewhere in the

34

 

 

 

 

material). In the pinning zone, the system is in a critical state, with supercurrent
density at the critical value. The pinning zone mitigates the current enhancing
effect of a crack by producing an effective crack tip radius of curvature of size p.
Provided p << a, we can combine the concept of a pinning zone with the explicit
calculations of Secs. [2.2] and [2.3] to predict the effect of a crack on critical
current. Using Eqs. (2.19) and (2.34) of Sec. [2.2] we thus find (for p < a < L
and ignoring logarithmic corrections)

jc(a) ~ 61(p/a)1/2, if ,0 < a < A (2 38)
jc(0) c2(p/A)1/2, if a >> A > p. '

where c; and c2 are undetermined constants that are independent of a. As is seen
from Eqn. (2.38), the main trend for jc(a)/jc(0) is that it reduces as the inverse
of the square root of crack size or penetration depth, whichever is smaller. This
prediction breaks down when pinning becomes sufficiently strong that the flux
pinning zone p becomes of order a, or L — a, whichever is smaller. When p is of
order L the critical state region extends across the sample, and we return to the

critical state model prediction of Eqn. (2.37).

The above predictions may be tested experimentally using superconducting
films. We suggest an experiment where a controlled crack of length a is introduced
into a high quality thin film (one in which we can control the pinning strength
relatively well, and which contains no extended defects which might compete with
the controlled crack). The crack can be drawn lithographically and the transport
critical current should be measured as a function of the crack length. We predict
that Eqn. (2.37 ) should apply if the film has sufficiently strong pinning (e.g.
highly polycrystalline Nb), while equation (2.38) should apply if the pinning is

weak (e.g. amorphous NbaGe).

 

 

 

Before closing this section, we make some comments concerning the effect
of finite film thickness and the attendant z-dependent current intensity on the
prediction (2.38). It is known that in films of finite thickness, there is both a z-
dependence and an z-dependence in the current density [Tink65, Duze81, Till86,
Orla91]. We are interested in the ratio jc(a)/jc(0), and we thus ask whether these
effects make a significant modification to the simplest estimates of this ratio.
Firstly the z-dependence. This is due to the complex form of the magnetic field
lines around the thin film, and may be calculated either within the Bean model
or within London theory. This effect will only modify (2.38) if the field profile
around the film is significantly different for a film containing a crack than it is for
a film without a crack. If the crack is much smaller than the film width (this is
the limit that we are primarily interested in), then this external field correction
is negligible, as can be qualitatively seen by considering the dependence of the
demagnetisation factor on crack length. This implies that the form of the z-
dependence of the current density is, at most, weakly dependent on the addition
of the crack. Assuming that the dominant part of the supercurrent is in the z-y
plane (as is expected for a thin film with current injected in the y direction and
the thin axis in the z-direction), we can use the theory of Secs. [2.2] and [2.3]
to estimate the current enhancement due to the crack, albeit with a 2 dependent
boundary field. Now, note that since the z-dependence is, to first order, only in ’

jg, the z-dependence drops out in ratio (2.18), and hence in (2.38).

A second possible complication for thin films is that there is an z-dependence
in the current density which is not given by the simple London form, but rather is
well approximated by j(z) ~ j(0)[1 — (z/L)2]1/2, for a: >> A and by an exponential
form for :1: < A. Since L is large compared to the a range (of order a) that we

36

 

 

 

are considering, we can consider the current flow to be essentially uniform for
the purposes of our calculations. This implies again, that the result (2.38 should

. provide a useful first estimate for a thin film, but now with a large effective A.

2.5 Conclusions

We have argued that there is a strong enhancement of supercurrent near cracks
which have a width greater than 6,, and length a >> width. We have illustrated
the effect using 2D London theory where a crack of length, a, leads to a strong
current enhancement as given in Eqs. (2.18),(2.19), (2.33) and (2.34). The main
scaling behavior contained in these equations is that the supercurrent increases
as the square root of the crack length, a, or the penetration depth, A, whichever
is smaller. Cracks can thus lead to large local current enhancements, especially

in geometries where the effective A is large.

We have also discussed the effect of large local supercurrents on critical current
measurements. We argued that large local supercurrents act as nucleation sites
for vortex penetration, and that these vortices feel a Lorentz force that tends
to make them move and hence cause dissipation. Balancing this, is the pinning
strength of the material. If the pinning is sufliciently weak, dissipation (vortex
nucleation and flow) is initiated by the large supercurrent at the crack tip. In
this case equation (2.38) should apply, and the superconducting material, is very
crack sensitive. If pinning is sufficiently strong, the vortices nucleated by the large
supercurrent at the crack tip are pinned and resist the strong local Lorentz force
at the crack tip. In this case the critical state model applies [see e.g. equation

(2.37)], and the superconductor is relatively crack insensitive.

37

 

 

 

 

To test these predictions, we propose that experiments which measure the
critical current as a function of the size of a controlled crack of length, a, and
width, 11) > 6., should provide a new and useful measure of the current carrying
capacity and reliability of superconductors. The thin film geometry (thickness
5 A) may be most useful, as there, current flows almost uniformly throughout the
cross-section of the film. A crack perpendicular to the direction of current flow
will cause a large amount of supercurrent to be channeled around the crack tip,
leading to strong enhancement effects as predicted by equations (2.18) and (2.34).
In addition, in thin films, the crack may be drawn in a controlled manner using
lithographic techniques. For weak pinning films, the simplest London theory [see
Eqn. (2.38)] predicts that the critical current reduces as the square root of the
crack length or penetration depth, whichever is smaller, provided a < L. For
experiments over a broader range of crack lengths, the square root law of Eqn.
(2.38) may need to be replaced by a more general exponent (i.e. 2: as defined
in the abstract not equal to 1 / 2). Experiments comparing crack effects in weak
pinning films (e.g. NbaGe), where the hotspot theory [e.g. Eqn. (2.38)] should
apply; with strong pinning films (e.g. highly polycrystalline Nb) where the critical
state model [e.g. Eqn. (2.37)] should apply, would be especially useful.

Finally, we note that the ideas of this study have many analogies in mechan-
ical fracture. Vortices are analogous to dislocations, and the pinning zone near
the crack tip is analogous to the process zone near the tip of a mechanical crack.
Controlled notch testing is standard in fracture and is analogous to the controlled
crack test proposed above. The analogy between mechanical fracture and critical
current in superconductors becomes closer as supercurrents become more inhomo-

geneous, and the pinning becomes weaker. In this limit, the current flow patterns

38

 

 

 

and critical current develop an extreme defect sensitivity that is characteristic of

fracture and other breakdown problems.

39

 

 

 

 

 

 

Chapter 3

Elasticity and fracture of
disordered networks

3.1 Introduction

Many alloys and composites are disordered. In fact no material is perfect and
the strength of a given material is very much dependent on these minor flaws or
defects. Thus understanding the effects of random defects and their statistics is
very crucial in analyzing the material strength. The effective elastic properties of
these materials can be calculated using effective medium theory and other “ho-
mogenization” methods. Random networks with “realistic” interatomic and many
body potentials are now routinely simulated on a computer. In combination with
perturbation theory near the pure limit, rigorous bounds, and scaling analysis in
the case of fractal geometries, there is now a well defined set of tools to accu-
rately predict the elastic properties of many random materials. With this success
there has been a growing desire to extend our understanding to the effect of ran-
dom microgeometries on the non-linear elastic response and especially the failure
strength of alloys and composites. It is known from the classic work of Griffith,

that penny shaped voids can act as nucleation sites for fracture. However, the way

40

 

 

 

in which crack nucleation and growth is affected by a disordered microstructure

is still poorly understood.

During the last few years, researchers, predominantly from the statistical me-
chanics and solid state physics communities, have introduced several new ideas
into the study of fracture. Firstly, the idea of a fractal geometry was applied to
fracture. surfaces[Mand84, Unde86, Mech89]; secondly fracture patterns were stud-
ied using the ideas of non-equilibrium growth [Tak386, Term86, Loui87, Meak88,
Hinr89]; and thirdly the methods of disordered systems were used to calculate the

average strength and statistics of random materials.

Due to the fact that fracture is initiated in regions of high stress, it is to be
expected that homogenization methods, unless they focus in the crack growth
region, should be poor predictors of the effect of disorder on failure strength.
Recently it has been shown that in random networks, failure strength is singu-
lar in the pure limit[Duxb87a, Duxb91, Bea188], so perturbation theory, in the
conventional sense, is unable to probe strength behavior. Useful bounds on ten-
sile strength have not been developed, probably because standard bounds lead
to strength limits that are so far apart that they offer no useful information.
Understanding of failure is then based on Griffith’s “single crack” ideas, ex-
tended in myriad ways, attempts to homogenize the microstructure (probably
incorrect in most systems) and more recently, detailed analytic and numerical
analysis of idealised disordered microstructures[Herr90, Sahi9l]. Along with new
scaling theories specifically taking into account the extreme statistics of brittle

fracture[Duxb87a, Duxb91, Bea188], these simulations have lead to new insights

into the fracture of materials.

 

 

 

 

 

 

 

In Sec. [3.2] the central force networks are discussed. The networks consist of
nodes and bonds. The nodes are connected by the simple harmonic springs (or
Hooke’s springs). The theory of elasticity based on the effective medium theory is
reviewed in Sec. [3.2.1] and the numerical simulation results using the conjugate-
gradient method is presented in Sec. [3.2.2]. In Sec. [3.3] the scaling behaviour

of the tensile fracture stress is developed based upon these results.

For some materials, such as graphene layers, the central force network can not
model the mechanical behaviour of the materials properly. In Sec. [3.4] I con-
sider the bond bending networks and their mechanical response by the numerical

simulation on honeycomb networks.
In Sec. [3.5] I summarize this chapter and draw conclusions.

This chapter contains materials which have appeared in the following three

publications:

0 [Duxb90a] P.M. Duxbury and S.G. Kim, Scaling theory and simulations of
fracture in disordered media, Proc. ASME (1990).

e [Duxb91] P.M. Duxbury and S.G. Kim, Scaling theory of elasticity and frac-
ture in disordered networks, Mat. Res. Symp. Proc. 207, 179 (1991).

e [Duxb94b] P.M. Duxbury, S.G. Kim and P.L. Leath, Size efl'ect and statistics

of fincture in random materials, Mat. Sci. and Eng. A176, 25 (1994).

3.2 Central force networks

Consider a regular lattice containing N central force springs which connect the

nearest neighbor sites of the lattice. Randomly remove volume fraction, f, of the

42

 

it

 

 

 

springs. As f increases, the elastic moduli decrease, upto the rigidity percolation
point[Thor87], f. (z 0.33 for the triangular lattice), at which the elastic moduli
vanish. If the springs are allowed to break when strained beyond a threshold
value, cc, the network maybe used to study damage evolution. In this case we
slowly increase the strain, allowing local spring rupture until the damage spreads
across the network and catastrophic failure occurs. The total energy, E, of a

spring configuration is given by

E = Z: ékijaij “10)2 (3'1)

where leg,- is the normalized Young’s modulus (spring constant) between adjacent
sites of the lattice, 1,3 is the length of the spring connecting the sites i and j, and
lo is the natural length of the springs (assumed constant). Disorder is introduced

through 13,-, which in the case of dilution are either 0 or he.

3.2.1 Effective medium theory for elasticity

The simplest estimate of p. = 1 — f. is found using a constraint counting
argument[Feng87, Thor83]. When the bond occupation probability, p, is small
the system consists of disconnected pieces and hence has many zero frequency
modes. The number of such modes being approximated by the number of degrees
of freedom (Nd) minus the number of constrains (zN p/ 2). Here, 11 is the spatial
dimension and z is the nearest neighbor co-ordination number of the lattice. Thus

the effective medium theory estimate of the fraction of zero frequency mode is

 

zsz—ZNp/2=1—f£

Nd 2d. (3.2)

43

 

 

 

When :1: —-> 0, the network is no longer “floppy”, and this defines the rigidity

threshold to be,
p. = 1 — f... z 2d/z (33)

Effective medium theory has been used to find effective elastic properties on ap-
proach to f... The spirit of the theory is to replace the random network with
a regular lattice having an effective spring constant km on each bond. This as-
sumes that the macroscopic symmetry is unaltered by disorder. Several different

methods have been used to estimate km and they all lead to the result

or: _ .
km '1 p p (3.4)

where ng are the elastic constants of the network, and the subscript or superscript,
m, denotes the effective medium value. Numerical simulations on triangular and
fcc lattice show that Eq. (3.4) provides an excellent approximation for the elastic
properties of random but rigid networks, for all but a very small region near the
rigidity percolation point (see also Fig. 3.6 for a comparison with the triangular

lattice result).

The tensile fracture stress of spring networks has also been recently studied,
and it is useful to first review the algorithm used in these calculations. It turns
out that there are many possible algorithms, and the choice of algorithm strongly
affects the final pattern of fractures springs. In contrast, the fracture strength

appears to be robust, and this is due to the fact that these networks are brittle.

 

 

 

 

3.2.2 Numerical simulation of damage and fracture

We apply an external strain to the networks and relax them using a conjugate
gradient method, with a convergence criterion of 10‘8 in the energy. When one
of the springs is strained beyond a critical strain, 6,, == (1;,- — lo)/lo = 10‘“, the
spring breaks (we set its kg,- to zero). In the results to. be described below, we
take a diluted triangular spring network, and slowly increase the external tensile
strain. At each external tensile strain we relax the network to minimize the total
energy of the system, Eq. (3.1) and check the strain in each bond against the
spring rupture criterion. If at a given external strain, a spring breaks, we fix
the external strain at that level, and continue relaxing the network and breaking
the spring carrying the largest strain (provided it is greater than cc) until either
catastrophic failure occurs, or until no more springs break. We call this the
hottest bond algorithm. We then slowly increase the external strain until it is
just large enough to break another spring in the network and again iterate the
hottest bond procedure described above. The stress—strain relationship of a single
spring, along with the stress-strain behavior found from these simulations on a
50x50 triangular network are shown in Fig. 3.1 and 3.2. It is seen from Fig. 3.2
that a catastrOphic event usually occurs very soon after the first spring in the
network breaks. We call the stress at which this event occurs o'1( f) A second
important stress is the stress required to rupture the entire network. This is the
maximum stress in the stress strain curve, and we label this 05( f) In networks
with small dilution (f small), o1(f) 5: 05( f) while at larger dilution this is not
always the case. Here, as in previous works[Bea188, Ha8589], we concentrate on

0'1( f), as this models the brittle fracture of the networks. The large strain behavior

 

 

 

-O

 

10"

Figure 3.1: The stress strain behavior of a single spring.

46

 

 

 

 

 

 

 

 

 

d
I ... 1.0-10'5
d
«1
C d
- 1040-5
4 .
d
. i
.. 4
o m :11 1010*
C ‘4
J
J
- 1.0104
d
(
ll/ ‘
#1; A L l 4L A A l L L 4L L .1 _jn L Hm Elm ‘) 0
0.0000 0.00“ 0.0000 Gm 0.01”
s

Figure 3.2: The stress strain behavior of 50x50 triangular network with f=0.1

47

 

 

 

 

involves a great deal of bond rotation, and is probably not relevant to either brittle
or ductile fracture. To illustrate the model network response, we show in Fig. 3.3
to 3.5 three snapshots of the microstructural evolution of damage and fracture
in a 50x50 lattice. The initial catastrophic event (usually at 01(f)), leads to a
path of broken bonds across the network, but not total fracture. Final fracture
occurs at much larger strain, though the final fracture path is often close to the
path of catastrophic event, and is usually preceded by bond rotation along the
fracture path. There are exceptions to this rule however, and crack branching
events do sometimes occur. As the initial f increase, the crack path becomes
more tortuous, though the structure of crack path is dependent on the model and
damage evolution algorithm used. The fracture strength however is more robust
and its scaling behavior appears to be fairly model independent[Bea188, HassSQ]
This makes the use of simple models most valid for the calculation of strength
and strength statistics, while the crack topologies found using these models are

only valid for the specific loading condition and microstructure that applies to the
model.

3.3 Scaling theory of tensile fracture stress

Since for small f, o1( f) e: ob(f) we can approximately calculate the fracture stress

from the relationship
660') N 00

05(0) ’ mm (3'5)

where 00 is the applied stress. omu( f) is the stress in the spring with the largest
strain, prior to any bond fracture. Eq. (3.5) uses the (small strain) linearity of

the system to extrapolate from a linear elasticity calculation of om( f) to deduce

48

 

 

 

 

 

 

 

 

 

 

 

Figure 3.3: The the network configuration after the first catastrophic event at
01(f). The simulation were on 50x50 lattices with f (initially) = 0.10

49

 

 

 

 

 

Figure 3.4: An intermediate configuration just before the second major peak in

 

 

 

 

 

 

 

 

 

Figure 3.5: The final fracture configuration. Full line from one node to its neigh-
bors represent the connected bonds while the truncated lines represent the bonds
which were initially present but which broke in the damage evolution process.

 

 

 

 

 

o7,( f). The problem now reduces to estimating om.x( f) for the initial disorder in

the network.

It is well known that cracks lead to large local stress enhancements. In dilute
networks, crack-like flaws are unlikely, but do occur with finite probability some-
where in the networks. It is these unlikely flaws that dominate the strength of the
networks. The simplest crack-like flaw contains n adjacent removed bonds. The

probability of occurrence of such a configuration is approadmately
P(n) c: Ldf" as n,L -+ co and f -—> 0 (3.6)

where L is the linear dimension of the networks in lattice units. The largest such

configuration that occurs with finite probability is estimated from

Ldfmw 2 1 (3.7)
which implies
nm 2 —dln L/ln f. (3-8)

Such a crack-like defect configuration has a stress intensity at its tips that scales
as (for —lnL/lnf large);

”Li—Q 2 1 + Km(-d1n L/ln n”? (3.9)
where K... is an unknown constant that depends on the curvature of the crack tip
and the size of a plastic zone (this is the continuum elasticity result for tensile

loading of a through crack).

In two dimensions (and for cracks for fixed finite separation), it is known that

two adjacent cracks enhance the stress between them by a factor approximately

.52

 

 

 

Ff

proportional to crack length[LiY887]. In this case the exponent in (3.9) is changed

to 1 [Bea188, HassSQ]. From these arguments and (3.5) we then find;
WU) ~ 1
do _ 1+Km(—dlnL/lnf)5’
for —d1nL/lnf —-> co and f —2 0 (3.10)

 

where 1/2(d — 1) _<_ B S 1. Eq. (3.10) also applies to the tensile fracture stress of
three dimensional system, with only minor modifications of K m and B. Eq. (3.10)

contains two interesting affects:

1. At fixed f, a logarithmic size effect in fracture stress occurs.

2. the fracture stress is rapidly decreasing for small 1", and in fact is singular

asf—20.

Eq. (3.10) is only strictly valid for —dln L/ ln f —+ 00, which does not apply to
the lattice sizes available for simulations. Despite this, the general trends outlined
in points 1. and 2. above are present in the numerical data presented in Fig. 3.6
and 3.7. For comparison the Young’s modulus which in this limit is linear in
dilution and clearly non-singular. The difference in dilute limit scaling behavior

of the Young’s modulus and the tensile fracture stress is clearly evident in Fig. 3.6

 

with the theoretical prediction (3.10) suggests that the scaling behavior of small

system may be well represented by the form;

01(f)~ 1
00 “1+Bllnfl-B

 

 

(3.11)

where B and B are system size dependent parameters. A test of this form using

the tensile fracture data of Fig. 3.6 is presented in Fig. 3.7. It is seen from Fig.

3.7 that B = 1.7 :l: 0.1 gives an excellent fit to the data in the singular region
53

 

 

 

 

 

 

 

J
1
.1
J
J
-l
J
1
l

 

vumnlmmam|aaanlnmmnlLamnlm

 

 

 

 

Figure 3.6: The Young’s modulus, Y(D), and tensile fracture stress, 01(f)(<>) of
L = 60 central force networks. Each point is an average over 40 configurations.
The solid line is from Eq. (3.4) applied to the triangular lattice.

54

 

 

 

 

 

2.50+++vTfifiv fra‘fT.w_c

 

 

 

r 4
_ 00 1
2.2: °
f’ ‘1
r °° 3° 1
C o o I
2 00 )— o o _.
r- o e 4
A " 4
=3, )- 0 ° 1: ° .
S 1 vs L- 0" 0 ° .1
9:. l o o o ,
f * 0° 0 '1
P 4
1 50 r- 59 0‘9 A
. a can I
O
i 0° 0 ° ‘
1 25 )- o o J
E 1
1 0° 1 L l J A EL L l . 1 A A l A J M J
0.0 0.1 0.2 o 3 0 4
twanv'P

Figure 3.7: A test of the scaling of ob(f) of the form (3.11) for 01(f) for B = 1.0
(0) and B = 1.7 (D).

55

 

 

 

 

f < 0.10. We thus suggest that Eq. (3.11) provides a useful representation of the
strength of the dilute random networks for a wide range of system sizes and that
B only lies in the range 1/2(d - 1) S B S 1 in the infinite lattice limit. We expect
that in the dilute limit, models with bond bending forces should show a behavior
similar to the central force networks. This may be understood if we consider a
lattice constructed of blobs of central force elements. Since the lattice is rigid on
most length scales near pure limit, most of these central force blobs are rigid. We
can thus renormalize the problem on to an effective model with bond bending
forces. In contrast, the behavior of the bond bending and central force problems
are very different at higher f. The bond bending problems lose rigidity at fc, the
connectivity percolation point, while the central force networks lose rigidity at f.

(f. < fc) [Thor87].
In bond bending problems, 05( f) exhibits the critical scaling;

M g (f. __ fy, (3.12)

‘70
Several workers have placed bounds on 2:, though its precise value is as yet un-
known. A node links and blobs picture[Guyo87] which states that the problem

may be replaced by a network with effective lattice constant £ (the percolation

correlation length) implies the upper bound,

a,( f) < sou. — no“): (3.13)

where u is the critical exponent for £. If the effective lattice constant, 6, acts as
a rigid moment arm on a given bond, the fracture stress will be greatly reduced.

This leads to the lower bound;

0'6”) > 00(fc — fldu- (3.14)
56

 

 

 

 

We thus deduce that z in Eq. (3.12) lies in the range (d — 1)u < a: < du. We
do not at present have any bounds on 2 for the central force problem, but it is
evident from Fig. 3.6, that :1: is nearly 1 for f quite close to f... Some experiments
have been attempted to find 2:. Using a 2mm thick aluminum sheet with holes
drilled to form a triangular lattice, Sieradzki and Li [Sier86] found a: = 1.7 :1: 0.1.
This value lies within the bounds given above (with u = 4/3 in 2-D). It must
be remembered though that the experiment is not truly in 2-D, the aluminum is
ductile and the lattice size is small (20x20) so this agreement may be fortuitous.
A second interesting feature of the Sieradzki and Li paper (see Fig. 2 of [Sier86])
is that it shows some experimental evidence of the qualitative difference between
dilute limit scaling of the elastic modulus and that of the fracture stress, despite
the small system size. A second experiment to find a: was published by Benguigui
et al [Beng87] who randomly punched holes in a copper sheet. They used a square
lattice and found a: z 2.5 :l: 0.4. This also lies within the bounds given above but

is also subject to the caveats stated above.

3.4 Bond Bending N etWorks

3.4.1 Background

In carbon fibers, the basic building blocks of the fibers are planar networks
made up of connected benzene rings. These graphene layers are imperfect
[Ke1181, Dres88, Donn90], and contain many vacancy and vacancy cluster de-
fects. In this and following chapter we attempt to assess the effect of atomic scale
defects on the in plane strength of graphite sheets and tubes. The recent reports

of “buckytubes”, which appear to be composed of nanoscale graphite tubes, lends

57

 

 

further impetus to this effort [Iiji91]. The model we use only considers the elastic
energy of the graphite sheets, and ignores the re-hybridization and charge transfer
that occurs near defects. The purely elastic model provides a reference against

which later calculations, which include these effects, will be compared.

The model and algorithm we use is as follows. The lattice is assigned an elastic

energy given by the Kirkwood energy;

I
E = '2- Z 0000' — 19,-)2
u

1
+5 2: B,,k(co.90.-,-,. - 6039?“)2 (3.15)

ijlc
i, j label nodes of a honeycomb lattice, while 1,,- (lg) is the length (equilibrium
length) of the bond between i and j and 01,-,- is the central force spring constant
between nodes i and j. 95,-), is the angle subtended by nodes ijlc, and (6:31: is the
angular spring constant associated with deviations in this angle from its equilib-
rium value 991'? This potential is similar to the Keating potential, which has been
successfully applied to covalently bonded materials and illustrates the points we
make in this work[Keat66, Kirk39]. We choose the Kirkwood potential for our
initial study, as previous studies of the effect of bond defects on the elasticity of
honeycomb structures have used this potential, and provide a check on our anal-
ysis [Zab086]. The change to a Keating potential will only make small changes to

the numbers calculated here.

For pure graphite, the ratio B /a ~ 1/4, and this is the value that we use in
our simulations. We assume that bond rupture only occurs in tension, and that
the strain at which this bond rupture occurs is 20%, as found from calculation

of the strain at the peak of the force-displacement curve of the carbon-carbon

58

 

bond (other strong atomic bonds have a similar value for the strain at the peak
force). When a bond exceeds this tensile strain, it is removed from the network.
No deformation is allowed out of the plane of the graphite sheet (no buckling),
and in the fracture simulations, the tube maintains the same average diameter

(Poisson ratio = 0).

To simulate tensile fracture, we carry out the following algorithm [Arca85,
Sahi86, Duxb87a, Duxb91, Bea188]:

1. Increment the applied tensile strain by a small amount.

2. Search through the bonds in the network to see if any bond carries a tensile
strain which is larger than the critical value of 20%. If any bonds exceed
this value, remove the one carrying the largest strain. Then return to the

start of 2. and iterate until no further bond failure occurs.

3. If no bond carries a strain larger than the threshold value, return to 1.

Practically, we find that it is possible to predict the applied strain at which the
next bond will fail from the previous relaxation, so the applied strain increments
are variable. This leads to a large increase in the efliciency of the algorithm. The
fracture algorithm is based on our ability to find the equilibrium strain configura—
tion for a honeycomb lattice described by the strain energy (3.15), with arbitrary
defect configurations. To find this lowest energy state, we use a non-linear conju-
gate gradient technique. This is composed of successive linear conjugate gradient
steps, with non-linear interpolation after each linear conjugate gradient step. The
relaxation is stopped when the largest residual force in the lattice is 10", which

corresponds to a residual energy which is always less than 10‘8. Relaxation of

59

 

ll

 

 

a 400 node lattice typically takes 1 minute of CPU time on a Convex CX240

processor.

The fracture algorithm described above is “quasi-static” because after each
bond failure, the network elastic energy is relaxed completely. The corresponding
experiment would have a very low applied strain rate. Use of a network and the
lattice energy (3.15) restricts the calculation to brittle materials, as no atomic
sliding is allowed in the model (no plasticity or dislocation creation). This is
complimentary to the molecular dynamics (MD) technique which corresponds to
very high strain rates (a typical MD simulation runs for less than a nanosecond),

but which does allow bond sliding and the plasticity associated with it.
3.4.2 Simulations of honeycomb networks

The elastic properties of the model represented by Eq. (3.15) on honeycomb
lattices, have been previously studied in some detail [Zabo86], and it is known that
the behavior near the percolation point is given by CH ~ )1. ~ (fc — f )T, where p
is the shear modulus and T = 3.96 £0.04. For the bond diluted honeycomb lattice
pc = 1—2 sin(1r/ 18) ~ 0.6527 (f6 = l—pc). Most simulations have been performed
for small ratios of B / a which are convenient to emphasize the difference between
the behavior of rigid lattices (e.g. triangular) and non-rigid lattices (e.g. square
and honeycomb). For B / a = 1/4 as is typical for graphite and related materials,
the dependence of the elastic moduli on random bond vacancies is presented in

Fig. 3.8

Using the failure algorithm, we find the stress-strain curve for a sheet with 5%
random bond defects to be as shown in Fig. 3.9. 06 is the tensile fracture stress of

this random honeycomb. The stress-strain curve of this disordered tube is close

60

 

 

 

 

 

 

inf a l f l f 1 f -)
r 1
0s)- 3 e
O o.)— 4— : p
:1 '1-
\ 1\
:3. o —o .9
0.4 l— . -' an
L o J o
s
01*— —1
. 3 .
o, . 1 . 1 1 1 1
'0.0 01 0.: 0.: 00

Figure 3.8: Cu and p for a 400 node tube as a function of the volume fraction of
removed bonds. 0?, and p0 are the values in the pure limit.

61

 

 

 

T
.l
J
1
1
.4
.1
._J
W
.4
-1
__1
.1
I

  

LllJlllllLllJlLllllllLLLl

 

TUlllrrlliirl11TrI1IITI'I1j

 

 

 

l

a
l"
)-
-)-—
)-
AL

00
0
.°
he
0
N
O
(D

m
12'

Figure 3.9: The stress-strain curve for a 900 node tube with 5% random voids.

62

 

1’

 

 

 

to linear elastic up to 03,. In Fig. 3.10 and 3.11 the initial defect configuration
used to produce Fig. 3.9, and the final fracture configuration produced by the

simulation are presented.

The simulations are carried out with periodic boundary conditions, so the
simulations are for tubes. Since f is small, the fracture path is almost linear,
as there is little crack deflection to take advantage of pre-existing flaws. Despite
this fact, the fracture stress is reduced appreciably from its pure lattice value.
From Fig. 3.9, it is seen, for example, that the fracture stress, 0’, is 2.52, as
compared to the pure lattice value of about 5.924. It is known that bulk graphite
. fails at a strain 6,, ~ 0.1 — 0.3%, while carbon fibers have failure strains that vary
from 1% (high modulus graphite fibers) to 5% (low modulus carbon fibers that
are composed predominantly of graphene sheets)[Kell81, Dre588, Donn90]. As we
shall discuss below, random defects cause 6., and ac to decrease as the lattice size

(sample size) becomes larger.

To illustrate the effect of atomic defect clusters on fracture, we show in Fig.
3.12, the response of a graphite sheet at 20% external strain, to a crack like defect
cluster containing 10 missing bonds. The stress enhancement in the bonds at the
ends of the defect cluster is 2.56. The applied strain required to break this flawed

tube is then reduced from its pure lattice value by this factor.

When bond defects randomly occur in a graphite sheet, the strength is reduced
by infrequently occurring large defect clusters. To calculate the average value of
this reduction for our honeycomb lattices, we average the value of 06 found from
curves such as Fig. 3.9 over many configurations. Results for a, for various values

of N and f are presented in Fig. 3.13.

63

 

 

 

In...

383?:3':§'§t§:-:-$:§:§§f§:

3.3%: 83?: garage:
,l;"-8.o. 33:3. 33%?
3§3§35. 5353:3 5353535353535353

35§5~.53533 35353035. 353

l y.

Figure 3.10: The initial configuration of tensile fracture of a hexagonal lattice
containing 5% void volume fraction.

64

 

 

 

 

)0...

,.-~...
':-::§:"’-i§i§§:3‘ '33

......
00000

5 = =. . z ' .
I: .3. E 0.0 Q ~ ..:.:.....
1° w

Figure 3.11: The final path of tensile fracture of a hexagonal lattice. The dotted
lines indicate the broken bonds.

65

 

UYY

    
  

. 10”

Figure 3.12: The response of a graphite sheet containing a crack-like defect cluster,
at applied strain of 20%

66

 

3.;
b
.1
q
-.
q
_
.1
d
-.
q
—
q
.1
q
.(
_
..
.1
.
..

 

 

 

3
t ‘ :
0.0 '— —:
o ?‘x d
O .0 x i
b 0.5 a, _
\ 3 03 x l
o : ° .
b 0 4 :- D x "i
: 8 3
0.2 E- : ._:
: I xxxxx :

0.0 M L I I l I I I I I I l I I I I I

0.0 0 1 0.2 0.3 0.4

Figure 3.13: The average fracture stress as a function of void volume fraction, f , for
networks containing 100(x), 400(0) and 900(0) nodes. Each point is an average
over 20 configurations at fixed f.

67

It is seen from Fig. 3.13, that there is a rapid drop in tensile strength at small
f for all values of N. A scaling analysis that we have previously applied to other
failure problems in both electrical and mechanical systems may be used for this
problem [Duxb87a, Duxb91, Bea188]. We find that as N —> 00, with f finite but
small,

0’c(f,N) 1

«72 ~ 1+ k(—1nN/1nf)a (3.16)

where k is a constant, typically in the range 0.1 - 0.6, while a is an exponent that
lies in the range 1/2(d — 1) < a < 1. The analytic analysis [Duxb87a, Duxb91,.
Bea188] shows that 0,; is logarithmically singular as f -> 0, although the system
size must be large for the analysis to apply. Another very important feature of
both the numerical and analytic analysis is that 0., decreases logarithmically with
N. In fact as N -> oo, 0': = 0. However, because the decrease is only logarithmic,
even with N = Avogadro’s number, the reduction in a". from its pure lattice
value, is only about a factor of 100. Nevertheless this factor of about 100 is
vitally important in the consideration of materials for structural applications. It
is also important to note that from a theoretical vieWpoint, the thermodynamic
limit (N = 00) is meaningless in brittle failure problems as there a}, = 0. In the
theoretical analysis N must always be finite, and often the interesting analysis is

the behavior at large N. Note that in the limit f —-> 0, (and large N), 011 is linear
in f (and obviously non-singular).

Fracture behavior as a function of flaw volume fraction is thus markedly differ-
ent than that of elastic moduli, and this difference has been related to the different
scaling behavior of the moments of the bond stress distribution[LiY889].

It is difficult to test the theoretically predicted logarithmic size dependence
68

 

 

 

in these graphite sheets as we are only able to simulate samples over two orders
of magnitude in N. Work on simpler models however allows a more detailed
theoretical and numerical analysis [Duxb87a, Duxb91, Bea188, Phoe92], which
confirms the logarithmic size effect, and suggests that it should be present in
models like that treated in our study. The analytic analysis also shows that
the variability in the mean of 0;, should be anomalous. In contrast, the elastic
moduli, for example the Young’s modulus, for two configurations at the same 1‘ of
a sheet with N nodes have values which are nearly the same: E1(N) ~ E2(N) +
0(1/\/N). However the same two configurations can have 0., values which show
much more scatter, and typically, 02(N) ~ 02(N) + 0(1/ In N). Thus although
06 self-averages to a unique value as N -> 00 it does so only logarithmically. In
contrast, the elastic modulus (away from pc) self-averages like a Gaussian random

variable.

3.5 Conclusions

A comparison of the elasticity and fracture of random networks illustrates the
different scaling behavior of these properties in the presence of disorder. Eqs.
(3.3) and (3.10) along with Fig. 3.6, show that the elastic constants are non-
singular except near special points such as the rigidity or connectivity percolation

points. In contrast, the brittle fracture stress exhibits the following features.

1. There is a dilute limit singularity in tensile strength (see Fig. 3.6, which
is well represented by Eq. (3.11). A dilute limit singularity is present in
some published experimental works (see e.g. Fig. 2 of [Sier86]), though
the specific form (3.11) has not been tested experimentally. It should be

69

 

 

 

 

emphasized that experiments on fracture strength trends must be performed
on several samples, as sample to sample variations can mask the predicted
trends in average strength. This is illustrated for example in the analogous
problem of dielectric breakdown in metal loaded dielectrics where the trend
analogous to Eq. (3.11) was observed, but about 10 samples at each value of
f were required to reliably see the trendICopp89, Duxb90b]. A logarithmic

size effect also occurs in the dilute limit [see Eq. (3.10].

. The statistics of fracture in the dilute limit obey a modified Gumbel form
rather than the usual Weibull distribution. Although many hundreds of
samples are required to differentiate between Weibull and modified Gumbel
statistics [Duxb87a, Bea188, Hass89], extrapolations to a reliabilities of order

10'6 leads to design predictions that can differ by order of 30%.

70

 

 

 

 

Chapter 4

.Melting of carbon fullerenes

4.1 Introduction

Recent success in synthesizing Cso “buckyball” clusters in macroscopic quantities
by Kratschmer et a1. [Krat90] triggered enormous interest of physicists, chemists
and materials scientists alike in search of new man-made materials with unusual
properties. The uncommon “hollow soccer ball” (or “fullerene”) structure was
originally postulated for carbon cluster by Kroto et a1. [Krot85]. For carbon alone,
a whole spectrum of structures has been proposed (An extensive survey of carbon
cluster literature up to 1989 has been given by Welter and Van[Welt89]), ranging
from giant hollow fullerenes (A good review of fullerenes has been presented in Ref.
[Cur191]) to elongated “bucky tube” or “nanotube” [Iiji91] and three-dimensional
graphitic membranes with “plumber’s nightmare” structures [Len092, Vand92].
Some of these structures are illustrated in Fig. 4.1. The icosahedral Can, called
“buckminsterfullerenes” or “buckyball”, is the most spherical molecule in nature.
Each of the 60 carbon atoms that are located at the vertices of the truncated
icosahedron are equidistant from the center-of-mass of the “buckyball”. The 0120

isomer in Fig. 4.1 can also be considered to be a short carbon nanotube: extension

71

 

[I

 

 

 

 

 

 

             

   

q ,7 ' "i.
.z,.. a,

 

w ,...... at». a... .. .............
as. ,......” Mags
9. ,.. w ,,.. a

\e‘

 

a:
.x( M
”if

  
 
 

» \.\....? .. rm . .. .2 «a _. .
,. WW waxy sausaysmw.” MW ..,. mam”... ”mm“... mum”. ...............

e 4.1: Carbon fullerenes. Cw “buckyball”, Cm; and Cm) isomers.
72

F'

 

ll

 

along the long axis of symmetry would yield a carbon nanotube similar to those

actually produced and studied by transmission electron microscopy (TEM) [Iiji91].

Due to the high structural stability[Cur188] and strong bonds with light mass,
carbon fullerenes have many unique properties that make them technologically
important as well as scientifically fascinating. Despite extensive studies over past
few years, almost no information is available about their equilibrium phases at
temperatures close to and exceeding the melting point of graphite. We expect the
phase diagram of these strongly correlated structures to be very interesting and

to consist of several well identifiable “molten” phases.

A second reason to study the high-temperature behavior of fullerenes is
related to their formation. While bulk quantities of fullerenes can now be
routinely synthesized in a carbon arc [Krat90], the microscopic mechanism of
their aggregation from gas phase is still the subject of a significant controversy
[Waka92, Waka93, Ebbe92, Kern92, Held93, Hunt94]. Since the fullerenes are pro-
duced in the high temperature processes, such as are burning or laser vaporization
of graphite, knowledge of the thermal stability will be useful for controlling the
synthesis of fullerenes. We also believe that a detailed study of the annealing
process may shed new light on this subject, since intermediate structures, which
occur during thermal quenching in the rare gas atmosphere, may also be observed

while heating up these structures to high temperatures.

A third motivating reason for this study are recent collision experiments which
indicate that fullerene molecules are extremely resilient and only fragment at en-
ergies exceeding z 200 eV [Camp93]. Since in the collision process the kinetic

energy is transferred into internal degrees of freedom as “heat”, additional infor-

73

 

 

 

mation about the intermediate structures of colliding fullerenes could be obtained

by investigating those of the superheated molecules.

In an attempt to elucidate the above three problem areas, we performed a
detailed molecular dynamics (MD) study of the melting and evaporation pro-
cess of three prototype fullerenes, namely Czo, 050, and 0240. In our simula-
tion, we investigated the response of a canonical ensemble of fullerenes to grad-
ually increasing heat bath temperatures using Nosé-Hoover molecular dynamics
[Nose84, Hoov85, Alle90], complementing recent microcanonical ensemble simu-
lations on Coo and C70 [KimE93]. Of course, the quality of simulation results
depends primarily on the adequacy of the total energy formalism applied to the
fullerenes. Since our simulation also addresses the fragmentation under extreme
conditions, the formalism used must accurately describe the relative stability of
fullerenes with a graphitic structure and their fragments, mostly short carbon
chains. The corresponding transition from .9122 to sp bonding can only be reliably
described by a Hamiltonian which explicitly addresses the hybridization between
02.9 and C212 orbitals. A precise evaluation of interatomic forces is required to
obtain the correct long-time behavior of the system. On the other hand, the long
simulation runs necessary to equilibrate the small systems during the annealing
process, and the large ensemble averages needed for a good statistics, virtually
exclude ab initio techniques as viable candidates for the calculation of atomic

forces.

For this reason, we base our force calculation on a Linear Combination of
Atomic Orbitals (LCAO) formalism which conveniently parametrizes ab initio
Local Density Approximation (LDA) [Hohe64, Kohn65] results for structures as

different as 02, carbon chains, graphite, and diamond [Toma91]. We use an adap-

74

 

 

Ill

 

 

 

 

tation of this formalism to very large systems, which is based on the fact that forces
depend most strongly on the local atomic environment and which has been imple-
mented using the recursion technique [Zhon93]. The force calculation can be per-
formed analytically to a large degree. This not only speeds up the computations
significantly, but also provides an excellent energy conservation AE / E f. 10'10
between time steps in microcanonical ensembles, and a linear scaling of the com-
puter time with the system size. Most importantly, forces on distant atoms can

be efficiently calculated on separate processors of a massively parallel computer.

In Sec. [4.2] I will introduce the tight-binding method as the simplified LCAO
method. I will summarize the LCAO formalism as well as the parameters used
for carbon.

In Sec. [4.3] I will review the recursion method as used for carbon clusters
and present the implementation to carbon tight-binding formalism. I will also
demonstrate the accuracy of this technique.

Sec. [4.4] is devoted to the molecular dynamics formalism. Nose-Hoover canon-
ical dynamics method will be briefly summarized.

Sec. [4.5] contains my main simulation results. Using the tight-binding
molecular dynamics with the recursion method, I simulate the melting of car-
bon fullerenes. I present snapshots of the intermediate phases and investigate
the fundamental driving force for the “melting” by analyzing the entropy. I will

identify a previously unknown “pretzel” phase of Cso.

In Sec. [4.6] I draw conclusions of the simulations and summarize.

This chapter contains material of the following two separately published

manuscripts:

75

 

 

ll

 

 

 

 

s [Kim894a] Seong Gon Kim and David Tomanek, Melting the fullerenes: A

molecular dynamics study, Phys. Rev. Lett. 72, 2418 (1994).

s [Kim894b] S.G. Kim, W. Zhong and D. Tomanek, Synthesizz'ng Carbon

F‘ullerenes, preprint (1994).

4.2 The tight-binding potential for carbon

The total energy functional for carbon which we used for the fullerene melting
process was proposed by Tomanek and Schluter [Toma91] and its implementation
using the recursion approach was fully discussed in Refs. [Zhon93] and [Kim394b].

It is briefly summarized in this and the following sections.

Recently, ab initio molecular dynamics simulations of liquid and amorphous
carbon with 54 atoms have been carried out using the Car-Parrinello method
[Gall89a, Ga1189b]. Nevertheless, first-principles studies are at present still lim-
ited by their heavy demand on computational effort, especially to do the detailed
molecular dynamics simulation in the thermal disintegration process of fullerenes.
The strong covalent nature and environment dependence of interactions between
carbon atoms make the classical potential approach very difficult. With two-body
and three-body interactions included, the Tersoff potential or similar approximate
schemes [Ters88, Ball90, Che191] have been tuned to reproduce graphite and di-
amond structures and dynamical properties very accurately. But such potentials
will have difficulties dealing with fullerene melting or collision process which in-
volves both high coordination number environments (graphite or diamond like)
as well as less coordinated structures like linear chains and the 02 molecule. The

tight-binding approach, which includes electronic structure explicitly in the Hamil-

76

 

 

 

Ill

 

 

tonian, has proven to be more successful in this respect[Toma91, XuCH92].

In LCAO methods[Root51, Root60, Mull49a, Mull49b], introduced to crystal
calculations by Bloch[Bloc29], the one-electron wave functions are expanded in a

basis of atomic orbitals {1%)}

I'M = go... l0...) (4.1)
with

H I49...) = 6:: |¢a) - (4.2)

The one-electron Hamiltonian for the electron systems can then be written as

H = Z CacIana + Z tgajgciacj'g. (43)
to:

iaifi

Here cl“ is the operator to create an electron orbital Iota) at site i and tgajp is the

hopping integral describing the hopping of an electron from orbital Ma) at site i

to orbital I053) at site j.

Slater and Koster have shown that the hopping integrals between two given
orbitals can be simply represented as a combination of h0pping integrals with c,
1r, and 5 symmetry, hence reducing the number of independent quantities [Slat54].

Moreover, it has proven useful to restrict intersite hopping to the nearest neigh-

 

bors only, in the tight-binding approximation. The tight-binding method is much
simpler than ab initio calculations, yet preserves the basic physics by allowing mix-
ing and intersite hopping of s and p electrons. Instead of computing the various

hopping integrals by first determining the true crystal potential using an ab initio

technique, we treat them as disposable constants chosen to fit accurate ab initio

77

 

 

 

 

 

calculations for selected structures at restricted symmetry points of the Brillouin
zone. In the following, we will use the tight-binding technique with Slater-Koster
parameters as an “intelligent interpolation” technique between results obtained
using ab initio techniques such as LDA. We adopt the hopping integrals from Ref.
[Toma91].

We have found it useful to separate the cohesive energy of the system (with

respect to isolated atoms) into two parts [Fou189], as
_ coh = Ebs + Erepo (4-4)

The first energy contribution, Eb” is the electronic band-structure energy, which is
nonlocal by nature and reflects the hybridization in the system. The second term,
Emp, contains the internuclear repulsion and all other corrections to Eb. such
as the closed-shell repulsion, exchange—correlation and energy double—counting
corrections. As shown previously [Toma91], the one—electron spectra of different
carbon structures can be obtained to a sufficient accuracy by mapping ab initio
LDA band structures for the corresponding systems onto a tight-binding Hamil-
tonian. In this procedure, the essential information about the electronic structure

and many-body effects in the system is kept intact.

The one—electron band—structure energy is given by

E... = :3 U: EN,-(E) as — gage. . (4.5)

Here, the summation extends over all atomic sites i, N,-(E) is the local electronic
density of states, and EF is the Fermi energy which is a global quantity. The

reference energy of an isolated atom is expressed in terms of the energy levels 6.,

78

 

 

 

 

 

and the corresponding occupation numbers 11,30 which satisfy the condition

Er
an = f... N,(E)dE. (4.6)

With this definition, Eb, is zero for both empty and full bands.

For carbon systems, we base our electronic structure calculation on a
Slater-Koster parametrized [Slat54] four-state (3, 1),, pg, 1),) tight—binding
Hamiltonian[Toma91]. The diagonal elements of the Hamiltonian are s— and
p-level energies E, and Ep. The off—diagonal matrix elements are the hopping

integrals with an exponential distance dependence
V0035) = VXSW‘U) (4-7)

where H labels the hopping parameter. S(r) is a scaling function of interatomic

distance r

5.0) = exp {up [1 — (3)] } exp {,... [(3% 453%} (4.8)

where r0 is the equilibrium nearest-neighbor distance in bulk diamond[Kitt86].
Unlike the previous scaling functions [Good89, XuCH92], n5 is not necessarily
same for different hopping parameters and we replaced the power-law dependence
of leading factor by the exponential to remove the divergence at small distances.
In our Hamiltonian, we consider those atoms as first nearest neighbors which are
closer than the cutoff distance (16 and we include interactions up to second nearest

neighbors. The parameters involved in the band structure energy are listed in

Table 4.1.

The second term in Eq. (4.4), Er», contains the internuclear repulsion and all

other corrections to E5. such as the closed-shell repulsion, exchange-correlation

79

 

 

 

 

1“

 

 

 

Table 4.1: Parameters for the band structure energy

 

‘
—7

 

 

 

 

1.546 A 2.16 A

E, E,
-3.65 eV 3.65 eV
v.2. v.2. v.2... v.2...
-3.63 eV 4.20 eV 5.38 eV -2.24 eV
nsw nspa' nppa nm
2 2 2 3
To Tc no
10

 

Li
-_7

 

80

 

 

 

and energy double—counting corrections. We approximate Erep by a sum of “quasi-

pairwise” functions Er,

I
Erep = Z: E,(r,-,-, 2,). (4.9)
"J

I used the word “quasi” because it is not really pairwise. There is a small contribu-
tion from coordination number [see Eq. (4.10)] which introduces a deviation from
a pairwise form. But in most cases, this deviation is negligible. In the summation,
each pair of nearest neighbor sites i and j is counted twice. The function E, can
be written as the sum of major contributions from the interatomic separation and

a small coordination number dependent contribution
E,(r.',-,z.-) = [E1(r,-,-) + E2(z.-)]¢(r.-,-). (4-10)
Here 03(1') is another distance scaling function with the form:
pm = exp [—m<1)'"° — (9%] (4.11)
To Tc
and z,- in Eq. (4.10) is the effective coordination number of atom 2' defined by

. _ ___1___. 4.12
z, _ g 1 + exp['r(r,-,-/r, — 1)] ( )

E1 is given by an exponential function of r and E2 is a. function of z,
E1(r) = VI exp(—r/r1)+ Vg exp(—r/r2) + V3exp(—-r/r3) (4.13)

E29) = E§(1 — 71-3). (4-14)
The parameters involved in the repulsive energy are listed in Table 4.2. This

repulsive energy will ensure correct total energies for Cg, carbon chain, graphite,

81

 

Ill

 

and diamond structures with a large range of bond lengths, if the band structure
energies are calculated exactly. It guarantees a smooth transition between different

structures and thus is suitable for molecular dynamics simulations.

82

 

 

 

Table 4.2: Parameters for the repulsive energy

 

 

 

 

 

 

171. me 1' 7
2 11 50 0.474
V1 V2 v, E3
10.34 eV 317.41 eV 5.0 x 109 eV -0.42 eV
7'1 7‘2 1‘3 Tr
0.9 A 0.35 A 0.05 A 1.95 A

83

 

 

4.3 The recursion method

The inclusion of electronic degrees of freedom allows for accurate and physically
more insightful calculations at the expense of a heavier computational load. Typ-
ically the most expensive part of the calculation is the diagonalization of the
total Hamiltonian matrix, which increases rapidly with increasing system size
N. The typical scaling is N3 or N2 log N, which makes calculation of a large
system very expensive. This computational requirement is not acceptable when
studying the dynamics of very large structures or performing molecular dynamics
simulations of the aggregation process. However, in a covalent system, a local
disturbance will only affect neighboring atoms, and its effect will decrease expo-
nentially with increasing distance. On the other hand, local properties are affected
mainly by the neighboring atoms. Most linearly scaling methods to compute the
electronic structure such as the recursion techniques[Hayd80] and the density ma-
trix approach[LiXP93], are based on this fact. In this Section we present a new
approach to calculating the electronic structure, which scales as N rather than as
N 3.

The most difficult part of the total energy evaluation is an efficient scheme to
determine the local density of states. While N;(E) in Eq. (4.5) should contain the
essential physics associated with the nonlocality of bonding, the exact function
is not very important since the band structure energy, given by Eq. (4.5), is an

integral quantity of N,(E).

As mentioned above, the idea underlying this approach is that the interactions
between two sites do not depend significantly on the bonding topology far away.

This approach has been used successfully in the calculation of the electronic struc-

84

 

 

 

 

ture of amorphous semiconductors [Joan74]. It has been proven [see, for example,
Eqs. (1.11) to (1.17) of Ref. [Hein80]] that the local density of states at the site

i = 0 can be given by
N0(E) = 111% (é) macaw +1.). (4.15)

The Goo element of the Green function matrix is given by the Dyson equation as

(pip)...

E - do —bl 0
—b1 E — a1 —b2

 

G00057)

-l

 

 

 

= 0 -b2 E-a2
00
= 1 b, . (4.16)
17—40“ 1 b2
E_a1—E—a:—...

The continued fraction coefficients an and b: are related to energy moments pa =
I“ dEE°No(E) of the local density of states at the site i = 0. These coefficients

are obtained by tridiagonalizing the Hamilton matrix of the system,

(4.17)

But if only the first few recursion coefficients are needed, it is more convenient
to use their recursion relations[Hein80, Kell80, Luch87]. The method involves
setting up a new orthonormal basis set la"), n == 0, 1,2, ..., the first of which lug),

we choose as the particular atomic orbital Ida) or linear combination of orbitals

85

 

 

 

for which we want the local density of states or the the atomic orbital of the atom

being considered. With given I110) and

|u_1) = 0
b3 = (uoluo) = 1 (4.18)
we generate for i = 0, 1, 2, . . .;
as = (nil Hrs Ins) ,
lui-H) = (HTS - 6;) lui) - 5? It‘s-1) ,
63.. = (14.114...) / 04-14) . (419)

Charge transfer between inequivalent sites in the structure is reduced by an on-
site Coulomb interaction which, in the mean-field approximation, can be mapped
onto a crystal potential. Variations of the crystal potential impose, to a leading
order, a rigid shift on the local densities of states. This shift of the core and
valence levels does not affect the crystal cohesion, yet is reflected in different core
level binding energies. In our calculation, we determine this shift by imposing a
local charge neutrality condition Eq. (4.6).

The main advantage of using the recursion method is that in most cases, we
only need the first few recursion coefficients to have a reasonable local density of
states from Eqs. (4.15) and (4.16). Furthermore we can increase our accuracy by

BimPly including more levels in the continued fraction.

Carbon is the one of the most computationally demanding elements in terms of
the number of recursion coefficients needed, due to the complexity of its bonding.
Directional bonding, i.e. the preference for definite bond angles for sp, sp2 or sp3

type bonding, requires that at least the four lowest moments [1,. be included in

86

 

 

 

 

 

the continued fraction. Furthermore, it is essential to treat the s, 11,, p, and p2
subbands in the density of states separately in order to represent correctly the

bond stretching and bond bending forces.

The fifth moment of the density of states depends only on the first and sec-
ond neighbors of any given site. The corresponding information, contained in the
recursion coefficients an, and bin: = 0, 1,2, is determined by tridiagonalizing the
corresponding small submatrix of the Hamilton matrix. Truncation of the contin-
ued fraction after b3 would lead to a set of 5—functions for the density of states,
corresponding to a small isolated cluster with many dangling bonds[Luch87]. A
physically more reasonable approach to describe very large structures is to embed
the small cluster in an averaged environment with similar bonding. This can be
achieved by attaching a Bethe lattice to the cluster [Joan74] which, to a large
degree, is equivalent to using

53
E — a” - t(E)’

 

1(a) = (4.20)

giving the square root terminator t(E) [Beer82] in the continued fraction in
Eq. (4.16). The terminating coefficients a” and b3 are easily determined in the
bulk limit from the lower and upper band edges. In our calculation, these val-
ues are taken from the linear chain, the graphite monolayer or bulk diamond,

depending on the local coordination of the site i = 0.

The essential information for the calculation of deformation energies is con-
tained in the recursion coefficients a,,, b3, with n S 17, 1] z 2. Attaching a square
root terminator to the continued fraction of Eq. (4.16) at this point would mean

to ignore the specific bonding topology and the corresponding connected loOps

Which modify the higher moments and continued fraction coefficients beyond 1).

87

 

 

 

We have found it useful to consider continued fraction coefficients beyond 1; = 2 in
the expression for the density of states. Rather than calculating these coefficients
explicitly by tridiagonalizing an accordingly much larger Hamilton submatrix, we
use pre-tabulated coefficients an, bf, (32,11,215) for the bulk structures as a “patch”
which we splice onto the continued fraction. The set of coefficients is taken from
the graphite structure. In order to avoid strong oscillations in the density of
states due to this “patch” and the square root terminator, we mix in the “patch”
coefficients with an, b: and the terminating coefficients as, bf, so that no strong
discontinuities occur in an and b3, at n = 1] and n = u [Luch87]. The continued

fraction coefficients in the “patch” 1) < n < V are given by

E: - 1[(1—' 13) +(1+sinzr£)a]
n —- 2 Sin 2 an 2 u

- 1 1m 1m
2 _ _ _ . _ 2 - _ 2
bn — 2 [(1 sin 2 )5» + (1+sm 2 )bu] (4.21)

with a = (2n - 1) - V)/(V — 17). While the above procedure gives an improved
description of the density of states, the higher “patched” coefficients do not affect
the lower moments of the density of states and hence have a small effect on the
band structure energy. It should be pointed out that our method of termination
is tolerant of a possible inaccuracy in the terminating coefficients and does not
induce unphysical oscillations in the local density of states. With this terminator,

the local density of states can be written

1 1

 

1.2 (4.22)
E " a0 " 1 bf

E—a1_E-02-tlEl

The validity of our formalism has been demonstrated by W. Zhong and et al.

 

[Zhon93],, as shown in Fig. 4.2. The densities of states for three most common

88

 

carbon structures, obtained using the recursion method, are compared with the
results of conventional tight-binding band-structure calculations which diagonalize
the Hamiltonian matrix directly. Electronic local density of states N (E) (solid
lines), integrated density of states I f; dEN (E) (dashed lines), and the band
structure energy —Eb,(Ep) as a function of band filling (dotted lines) is plotted for
different carbon structures. Results of a tight-binding band-structure calculation
for an infinite carbon chain (bondlength doc = 1.286 A) (a), a graphite monolayer
(doc = 1.418 A) (b) and bulk diamond (doc = 1.546 A) (c) are compared with
the results obtained with our simplified recursion method in (d), (e) and (f),
respectively. Even though the recursion method misses some of the details in
the local density of states, it is very evident that the integrated quantities, such
as the band structure energies, are reproduced within excellent agreement. The
interested reader should be referred to Ref. [Zhon93] for a more complete list of
examples which demonstrate the sound basis for using the recursion method for

carbon clusters.

4.4 Molecular dynamics

Once the total potential energy calculation scheme is constructed using the recur-
sion method, the next step is to find the equations of motion according to which

the system evolves.

The dynamics of an isolated (microcanonical) N ~particle system in three-

dimensional space is governed by the Lagrangian
N 1
L = 23 5mm? - V({q.-})- (423)
5:1
Here, q.- is the position vector of atom i in the system and V({q,-}) is the total po-
89

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1.25 :- 1.25
1.00 g) 1.00
I: 0.75} 0.75
.2. 0.50: 0.50
5’3 0.25: 0.25
5, 0.00: , 0.00
'5 1.25;- 1.25
; 1.00:- 1.00
g 0.75: 0.75
0.50; / ’, ,—--' 0.50

O. i — .
25= ...... V‘. .rij 3::

 

0.00 ' ‘ ‘ ‘ ‘ ‘
-20 -10 O 10 20-20 -10 0 10 -20 --10 0 10
E [eV] E [eV] E [eV]

Figure 4.2: Comparison of recursion method vs. conventional tight-binding
method. Electronic local density of states N (E) (solid lines), integrated density of
states If; dEN (E) (dashed lines), and the band structure energy -Eb.(EF) as a
function of band filling (dotted lines) for different carbon structures. Results of a
tight-binding band-structure calculation for an infinite carbon chain (bondlength
doc = 1.286 A) (a), a graphite monolayer (dcc = 1.418 A) (b) and bulk dia-
mond (doc = 1.546 A) (c) are compared with the simplified recursive results in
(d), (e) and (f). The energy zero coincides with the Fermi level for a half-filled
band, corresponding to neutral carbon. (Fig. 1 from Ref. [Zhon93] used with the
permission from the authors.)

90

 

 

 

[001/110] ""3- 'soa 17°1dean

 

 

tential energy of the system. In present case it is given by the tight-binding Hamil-
tonian. The dynamical evolution of the system is described by Euler-Lagrange
equations derived from the above Lagrangian. This procedure yields statistics for
a microcanonical ensemble. Since we are studying the melting process of carbon
fullerenes, we want to fix the “temperature” of the system to a given heat bath
temperature. We define the temperature in such a small system by its instanta-
neous kinetic energy. We need to adapt molecular dynamics algorithm so as to
sample a constant-temperature ensemble. Several different methods of prescrib-
ing the temperature in a molecular dynamics simulation exist. A recent review
[Ande84] has attempted to summarize these methods and highlight their advan-
tages and disadvantages. Here I will discuss only one of these which I used for

melting the fullerenes.

We used the extended system method, or so-called N osé—Hoover method. This
method treats the dynamics of a system in contact with a thermal reservoir by
including an extra degree of freedom which represents that reservoir. I carried
out a simulation of this ‘extended system’. Energy is allowed to flow dynamically
from the reservoir to the system and back; the reservoir has a certain ‘thermal
inertia’ associated with it, and the whole technique is rather like controlling the

volume of a sample by using a piston which has certain mass.

Nosé achieved a major advance in this method by showing that the canoni-
cal distribution can be generated with smooth, deterministic, and time-reversible
trajectories[Nose84]. To do this he introduced a time-scaling variable s, its con-

jugate momentum p,, and a parameter (Nosé—Hoover parameter) Q. Nosé wrote

 

 

 

 

 

 

 

an augmented Hamiltonian

HNosé = 2P2/2m32+V({qa})

+pf/2Q + (X +1)kT1n 3, (424)

where X is the total number of degrees of freedom of the system. This Hamil-
tonian contains a nonlinear collective potential in which the time-scale variable 5
oscillates. Thus the system is coupled to a heat bath (described by the variables
s and p.). N osé proved that the microcanonical distribution in the augmented set
of variables is equivalent to a canonical distribution of the variables q,p’, where
the p’ are the scaled momenta p/s. Thus the Hamiltonian (4.24) generates the

canonical probability distribution independent of the values chosen for HNose’ and
Q.

The equations of motion from Nosé Hamiltonian (4.24) are

d = P/ma2
:3 = F(q)
5 = p./Q
p‘. = 2p2/m23—(X+1)kT/p. (4.25)

Hoover made a great improvement in implementing Nosé’s method computa-
tionally, by deriving a simpler form of the above coupled first-order differen-
tial equations[Hoov85]. He noted that if the time scale is reduced by s, then
dtold = sane... All of the rates given in Eq. (4.25) can then be expressed as
derivatives with respect to tn" (for which we will still use the superior dot nota-

tion)

q = p/ms
92

 

 

 

 

ll

 

 

 

 

 

 

 

15 = sF(q)
5 = spa/Q
p'. = 2:112/ms2 — (X +1)kT. (4.26)

The somewhat inconvenient variable s can then be eliminated from Eq. (4.26) by

rewriting the coordinate time evolution equations in terms of q,é, and 5:
31' = fi/ma - (P/m3)5/3

= F/m—éps/Q

F(q)/m— 44 (4.27)

The thermodynamic friction coefficient C E p, / Q, which appears in the second-

order Eqs. (4.27), evolves in time according to a first-order equation

6: [Z mcjz _ (X + 1)l¢T] /Q (4.28)

To have set of first-order equations, we redefine p E mq’ and replace Nosé’s X + 1

beobtaining
ti = 17/70
13 = F(9)-(p
( = [sz/m—XkT] /Q. (4.29)

This is the so-called Nosé—Hoover equations of motion. Hoover proved that the dis-
tribution resulting from Eqs. (4.29) is canonical in the variables q and p [Hoov85].
It can be shown easily that we can arrive at the above equations by generalizing
the microcanonical Lagrangian (4.23) to the augmented Lagrangian

N 1

L = 23 51.4.4253 — V({q.-}) + $052 — (X + 1)Tln , (4.30)

i=1

93

 

 

 

 

 

 

and applying Euler-Lagrange equations:

.4 3L. _ ab _ 0
dt 89's}: aqm — ,

4626:) = ..

The parameter Q controls how fast the system explores the available phase

 

space. For a large Q the trajectories gradually fill the phase space as expected
for an ergodic system, and for a small Q the trajectories develop more singular
turning points in phase space. For very large Q these equations simply reproduce
the microcanonical behavior. The proper range of Q values should be determined

by actual simulations for the given system at hand.

4.5 Simulation and Results

In our molecular dynamics simulations, we use the leapfrog technique to integrate
the equations of motion over time steps At = 5 x 10"16 s. We found that the
temperature of the system is controlled reasonably well when we use the value for
the Nosé—Hoover parameter Q = 1 / 250 (amu-A). This value corresponds, in crude
sense, to using imaginary thermal particles which are almost 3,000 times lighter
than carbon atoms as the mediator or carrier of the heat between the system and

the heat reservoir.

We begin each MD run by equilibrating the particular fullerene for over sev-
eral thousand time steps at a temperature T,- = 200 K, after the directions of
the initial atomic velocities have been randomized according to the Maxwell-
Boltzmann velocity distribution at that temperature. We increase the temper-

ature of the heat bath in steps of AT = 400 K from T.- to the final temperature
94

 

 

 

 

T, = 10,000 K, and let the fullerenes equilibrate for 800 time steps at each new
temperature. Each simulation run consists of more. than 20,000 time steps and
takes several hours of CPU time on a high-speed single-processor workstation.
Statistical (time-averaged) data for the structure and energetics are collected af-
ter the system has adjusted to the new temperature, which occurs typically after
400 time steps following a temperature increase. This procedure is necessary to
remove any undesired oscillation in the data due to the sudden increase of the
heat bath temperature. A separate ensemble averaging at each temperature is
necessary since ergodicity is not guaranteed especially in very small systems. We
perform ensemble averaging over 50 complete runs with different initial states of
020, 50 runs for Ceo, and 15 runs for 0240. This reflects the decreasing fluctuations

of thermodynamic quantities with increasing system size.

In Fig. 4.3, we present results for the temperature dependence of the total
energy E per atom for the 060 molecule. We observe a steady increase of the total
energy E with increasing temperature, with a well-formed step at T m 4,000 K
indicative of a “phase transition”. To get a better understanding of the structural
transformations in this and other fullerenes, we investigated the specific heat
not only of the Can, but also the 020 and the 0240 molecules. The calculated
temperature dependence of the specific heat per atom of these systems, obtained
from cv .= dE/dT, is shown in Fig. 4.4.

In the temperature range addressed in this study, we found significant devia-
tions'from the classical value cv = 3kg to occur only for Th2,000 K. We found

all features in the rich structured cv(T) to be reproducible from run to run up

to a temperature of T m 6,000 K. Individual peaks in cv reflect the latent heat

of transition between different “phases”. As we discuss in the following, these

95

 

 

 

 

 

 

l
1!
C
d
.1
d
.(

—-1
.l
d

 

 

 

q
)- c1
1- ‘-

-4.5 '—
- q
S‘ E I
u- d
.2. -5~° .— ‘1
fl " «I
g : :
g -5 5 _— -
a 0 0 : :
- — —
19- ° - -
l- q
,. a:
-6.5 _ _
d
d
-7 o I I L I I I I I I I L I l I I I I I I I I I ‘

 

E

2000 4000 0000 0000
Heat Beth Temperature '1' [I]

Figure 4.3: Temperature dependence of the total energy E per atom for the 060
fullerene (solid line). The data points in the E(T) curve for the fullerene mark
the discrete steps of the heat bath temperature used, and the dashed lines with a
slope of 3103 are guides to the eye.

96

 

 

 

 

 

     

 

 

 

pi I—l fil rm I T’; III l—‘Wfi l I I I T T
I
5.0" / l‘
.- I
I \ I
4.01- ’ 1 ..
l
in ’~’ \ /\ l \
3.0 s o sssss ssssss sss 00’ssss‘ssssssssad ssssssss Y 0“
b \ \ I ~
0 up
*0
3 4.0
:13
g 5.0
8' 0.0r , 1.
als‘". I
I- C20 0‘. ,.‘. o’0 ‘.’s‘. ,..]
3.0pssssss sssss so sssssss s sssssss .qu"lungs—s3...”Yunnan........\..\.1
s I ° 0’ ‘0’
—.- .I \o-.ao‘ ’0 \.’ —-1
I .\ I .\ s
2.0" ' v .
1 LI I I I I I I I I_ I I j 14 Li I J L
2000 4000 6000 8000 10000

Heat Bath Temperature ‘1' [K]

Figure 4.4: Specific heat cv(T) per atom of icosahedral fullerenes Cm (dashed
line), Coo (solid line), and ON (dash-dotted line), in the temperature range T <
10, 000 K, in units of leg. The ordinate is displaced vertically between successive
plots for clarity; the classical value cv = 3163 is indicated by the dotted lines.

 

 

 

 

 

“phases” can be characterized by their topology and atomic binding energy dis-
tributions. In contrast to phase transitions in solids, the analogous transitions in

fullerenes are not sharp owing to the finite size of these systems.

Due to the similarity of the cv(T) results for the different fullerenes, we con-
centrate in the following on the 060 molecule. Each peak in the cv(T) curve
corresponds to the transition from one intermediate “phase” to next one. In order
to identify the characteristics of its individual “phases”, we first take the “snap

shots” of the 060 molecule at selected heat bath temperatures as shown in Fig. 4.5.

To have some quantitative criterion to characterize the “phase” we also plot the
atomic binding energy distribution in Fig. 4.6 at the temperatures chosen for Fig.
4.5. The atomic binding energy for a given atom is defined as the contribution of
the given atom to the total cohesive energy of the whole system. Thus the cohesive
energy in Eq. (4.4) can be written as

N
Ecoh = 2:1 Eb“) (4-32)

where the binding energy of atom i is given by

135(5) = [[l: EN,‘(E) dE - 2115,06,, + 2,13,.(1‘53', 2;), (4.33)

with E,(r,-,-,z,-) defined in Eq. (4.9). We find the binding energy distribution
to provide a better signature of the different “phases” than topological quantities
(such as coordination numbers) which depend on the definition of cutoff distances.
We also monitored the change in the distribution of bond lengths for each atom

and plotted them in Fig. 4.7 at the same temperatures as in the previous two

figures. Only the atom pairs within the cutoff distance (2.0 A in this case) were
98

 

 

 

 

(a) Solid phase (b) Floppy phase (c) Pretzel phase
T-1000 K T-SOOO K T-‘ZOO K

 

 

 

 

] (f) Chain gas
T-IOZOO K

14”.?“

T-5000 K

 

 

 

 

 

Figure 4.5: “Snap shots” illustrating the geometry of a 060 cluster at temperatures
corresponding to the different “phases” discussed in the text.

99

 

 

 

 

 

Probability Distribution

 

 

 

 

(e)

‘1'II4200K

(a)

“II-1000 K

 

 

 

 

  

  
  

  

 

   

(1)

1.10200

(d) (e)

r-sooox 5 r-ssoox

 

 

 

  

 

 

5 6 7 5’ 6 7
Atomic Binding Energy E00,, [eV]

Figure 4.6: Atomic binding energy distribution in a Cw cluster at temperatures
characteristic of the different “phases” discussed in the text and in Fig. 4.5. For
each structure, the average binding energy is given by the dotted line.

100

 

 

 

 

 

 

 

considered in Fig. 4.7. In Fig. 4.8, we plot the distribution of bond angles at
the selected temperatures. If an atom is connected to more than one neighbor
within the cutoff distance, the angles between these bonds were calculated and
their distribution is collected. As the final monitoring quantity, we calculated the
mass spectrum of the fragments of the C60 molecule at selected time steps. The
atoms are considered to be part of the same cluster when they are connected by
bonds within the cutoff distance (2 A). The changes of the mass spectrum as a
function of the heat bath temperature are depicted in Fig. 4.9. When the system
consist of one (connected) cluster, this spectrum will show only one peak at the
value of the total mass. When the structure is fragmented into smaller pieces,

several peaks will occur at lower mass values.

In the solid phase, occurring at T32, 400 K and depicted in Fig. 4.5(a), the 060
is intact. As shown in Fig. 4.6(a), the binding energy of all atoms is approximately
the same, but decreases gradually with increasing temperature due to thermal
expansion. The bond length distribution shows a very prominent single peak
around 1.46 A [Fig. 4.7(a)]. The experimental values of two bond lengths for 060
at ground state are 1.40 A (double bonds) and 1.45 A (single bonds) [Yann91].

A peak in the cv(T) curve at T as 2,400 K indicates the gradual onset of
the floppy phase, depicted in Fig. 4.5(b). While the system is topologically intact
in this phase, as documented by a relatively narrow bond length distribution
[Fig. 4.7(b)], we observe the bond angle distribution to spread out [Fig. 4.8(b)],
resulting in a significantly broadened binding energy distribution [Fig. 4.6(b)].
Assisted by local shear motion, carbon pentagons and hexagons tilt easily with
respect to the surface normal especially in large fullerenes with flat faces. In 060,

this floppy motion creates a distribution of carbon sites which are either closer to

101

 

 

 

 

 

 

 

(b) Floppy phase (c) Pretzel rinses
r-sooo K E r-ssoo I ' '

 

   

151000 I

(d) linked chains (0) Fragments
r-sooox : r-ssoox :

Probability Distribution

  

 

 

 

 

lat. II IILIIILIJL‘ II Ajl

0.5 1.0 1.5 0.5 1.0 1.5 0.5 1.0 1.5
Bond Length [A]

Figure 4.7: Bond length distribution in a Coo cluster at temperatures characteristic
0f the different “phases” discussed in the text and in Fig. 4.5. For each structure,
the average bond length is given by the dotted line.

 

102

 

 

 

nosgnsmfln haaadn°hm

 

Probability Distribution

 

(a)

1'10” K

 

 

(d)

 

(b)

(a)

m:

(c)

WK

(1)

1'10“ I

 

 

 

.M.

100 150

50 100150

We

50 100150

Bond Angles [Degree]

Figure 4.8: Bond angle distribution in a Cm cluster at temperatures characteristic
of the different “phases” discussed in the text and in Fig. 4.5.

103

 

.Iuln .Iqllu a
an.“
no:

01".“
can.
\Hfia

and

 

Probability Distribution

 

 

(a) 0:) (c) w
r-1000 x r-sooo x r-ssoo r
(a) (.) 3(1)

r-sooo x 5 r-uot x r-10s00 x

 

20 40 80 20 40 00 A 20 40 00
Cluster Mass [Carbon Mass Unit]

Figure 4.9: Mass spectrum in a 060 molecule at temperatures characteristic of the
different “phases” discussed in the text and in Fig. 4.5. The cutoff length 2.0 A
is used to determine the cluster boundary.

104

 

 

 

 

or further away from the optimum graphitic bonding geometry than the average
structure. We found that the local shear motion leads to substantial cluster shape

deformations.

Above T as 4, 000 K, we observe a dramatic transition to the pretzel phase,
consisting of interconnected carbon rings and depicted in Fig. 4.5(c). They are still
all one single connected body as depicted by a single peak in mass spectrum [Fig.
4.9(c)]. A similar, yet not as complex structure has recently been proposed as a
high temperature phase of On clusters based on diffusion data [Held93, Hunt94].
In our calculation, we find this structure to be initiated by the rupture of a sin-
gle bond connecting a pentagon to a hexagon, which creates a large Opening at
the surface. As temperature approaches to this transition temperature, this bond
opening is amplified and propagates in a “zipper” motion through the molecule.
Suddenly the closed structure “unwraps” to become a “pretzel” structure. As
shown in Fig. 4.3, the transition from a (two-dimensional) fullerene to a (one-
dimensional) linked chain structure is clearly reflected in a sharp increase of the
entropy towards a value characteristic of a (metastable) Cm chain. The ener-
getically less favorable sp bonding of a growing number of two-fold coordinated
atoms in the “pretzel” phase is reflected in the broadening of the binding energy
distribution especially towards lower values. This is shown in Fig. 4.6(c) for a
fully developed “pretzel” at T = 4,200 K. One thing to note about this distri-
bution is that it has two major peaks. The atoms which are part of the chain
is mainly two-fold coordinated and has smaller binding energies and responsible
for the peak at lower value. On the other hand, the atoms which are part of the
connecting jOints of the rings have three or even higher coordination numbers and

contribute to the peak at the higher binding energy value. The specific heat data

105

 

 

 

 

of Fig. 4.3 suggest that the critical temperature shifts upward with increasing
fullerene size. The bond angle distribution shows a lot of new peaks over 150

degrees representing the chain structure [Fig. 4.8(c)].

Closed carbon chains of the Geo “pretzel” Open up at T2.5,000 K, as shown
in Fig. 4.5(d). The number of structural constraints in this linked chain phase
decreases, allowing individual atoms to relax. This leads to a smaller number of
inequivalent sites and more pronounced peaks in the binding energy distribution,

depicted in Fig. 4.6(d).

Starting at T m 6,000 K, we observe fragmentation of all fullerenes which
we studied. The mass spectrum finally shows split peaks at the values other
than 60 [Fig. 4.9(c)]. Typical fragments, such as those shown in Fig. 4.5(e) for
T = 5, 400 K, have “pretzel” and “linked chain” structures discussed above. Above
this temperature, the thermodynamics of the system is that of fullerene fragments;
strong run-to-run fluctuations are caused by the differences in the specific heat of
the individual fragments. While the binding energy range is similar to that in the
“linked chain” phase, the distribution is smoother at this higher temperature, as
shown in Fig. 4.6(e). Finally, at temperatures close to T m 10, 000 K, a conversion
to a chain gas is completed, as illustrated in Fig. 4.5(f). The mass spectrum shows

many peaks at the small mass numbers [Fig. 4.9(f)].

The temperature scale relevant to structural transitions in fullerenes can be
linked to well established thermodynamic data for graphite [Weas90]. We find
the “floppy phase” of fullerenes to occur at temperatures close to the melting
point of graphite, Twp, = 3,823 K [Weas90]. On the other hand, the calculated

fragmentation temperature T as 5, 400 K lies close to the observed boiling point

106

 

 

 

of graphite, T1,,“ = 5, 100 K [WeasQO].

In order to understand the rich “phase diagram” of fullerenes, we have to in-
vestigate the free energy F as a function of temperature. It is intuitively clear
that the more “floppy” chain structures have a larger entropy and hence should
be preferred over the more compact fullerene structures at sufficiently high tem-
peratures. As mentioned above, we have determined the entropy per atom S (T)
from MD simulation runs of C50 fullerenes and metastable Cso chains. Our re-
sults, shown in Fig. 4.10, indicate that the entropy values are equal at very high
temperatures where also the fullerenes have converted to linear structures. Prior

to the conversion, we find for the entropy difference
S(Cao chain) - 5(060 fullerene) a: 2kg. (4.34)

This allows us immediately to estimate the conversion temperature from fullerenes

to chains from low temperature data. Noting at the moment of conversion

F(chain) — F(fullerene) = [E(chain) - E(fullerene)]
-—Tc[S(cha1n) - S(fullerene)]

= o, (4.35)

we have
T, = AE/AS. (4.36)

Using the result obtained previously[Toma91]

AE = E(chain)-E(fullerene)
z. 1 eV (4.37)

107

 

 

 

 

 

 

 

  
  
 
  
 

 

 

.l I I I l j T I I T I I l I— j j I l T T I l
E I
s L- ...
.. .l
2 - .l
01 '- .."'. a
E 4 L— ..." .-
a I chain 2
— fullerene -
2:- -i
. J
" 1
o #1 l I I I I I I I I I L I k I I I I I I L

 

2000 4000 0000 0000 10000
HutBath‘runper-atnrefm

Figure 4.10: The entropy 5 per atom for the C60 fullerene (solid line) and the
metastable Coo chain (dotted line).

108

 

 

 

 

 

g

 

 

 

and assuming that the values for AE and AS = S(chain) -— S(fullerene) are
nearly temperature independent, we obtain Tc x 5,800 K for the fullerene-to-
chain conversion. This is in remarkable agreement with our simulation results

considering the crude arguments we have used in this estimate.

It is interesting to note that the open fullerene phases, which were identified
in our simulation, also occur during energetic collisions between fullerenes. In
particular, the final state of a Coo cluster undergoing an inelastic collision with a
0240 cluster at 300 eV center-of-mass kinetic energy heats up during the collision
process to temperatures close to 5,000 K and shows a final structure similar to

that presented in Fig. 4.5(d) [KimS94b].

Our results shed new light also on the formation process of fullerenes from
the gas phase. While time reversal of the evaporation process studied here would
provide one possible aggregation path, the general formation mechanism is clearly
more complex. We expect that many structures which occur during the frag-
mentation also reoccur during the aggregation. Of course, most of the “random”
aggregation paths will not lead immediately to a closed fullerene, yet are likely to
contain stable structural elements such as chains and rings. Once such a “linked
chain” or “fragmented pretzel” structure is formed, the aggregation dynamics will
be governed by many unsuccessful attempts to reversibly “roll up” the chains to
a fullerene, followed by one successful attempt to create an inert structure. This
picture agrees with the fullerene formation mechanism suggested in Ref. [Hunt94]
and the isotope scrambling results of Ref. [Ebbe92] suggesting that gas phase

assembly of fullerenes starts from atoms and very small carbon clusters.

109

 

 

 

 

4.6 Conclusions

We presented the new and efficient recursion method applied to the tight-binding
molecular dynamics for carbon atoms. The ’carbon potential parameters fitted to
LDA calculations have been presented and the reliability of the recursion technique
treatment has been demonstrated. Finally, we presented molecular dynamics sim-
ulation results of the melting and evaporation of the carbon fullerenes Czo, 060,
and 0240. Among others, we found a most dramatic transition to a “pretzel” phase
to occur at a high temperature of TZAOOO K. The temperature where fullerenes
disintegrate into carbon chains was explained using quantitative results for the

entropy.

110

 

a?"

 

 

 

 

 

 

 

Chapter 5

Stability of Magnetic Dipole
Structures in Ferrofluids

5.1 Introduction

A ferrofluid is a magnetic colloid composed of magnetite particles suspended in a
liquid, usually a petroleum oil [Rose85, Wang94]. In many ways, ferrofluids be-
have similar to electrorheological (ER) fluids which consist of a suspension of fine
dielectric particles in a liquid of low dielectric constant [TaoR92, Bloc87, Fili88].
The aggregation of particles in ferrofluids and ER fluids has been the subjects
of many experiments and numerical simulations. It has been observed that the
aggregated Co particles, produced by inert-gas evaporation in Ar, show very dif-
ferent fractal dimensions depending on the strength of the magnetic interactions
between the particles[Nick88]. Remarkable labyrinthine patterns are formed when
a droplet of ferrofluid is trapped between two horizontal glass plates in a vertical
magnetic field [Lang92, Dick93]. The effective viscosity of an ER fluid increases
dramatically if an electric field is applied, and when the field exceeds a critical
value, the ER fluid turns into a solid. Experiments also find that upon applica-

tion of electric fields, dielectric particles in ER fluids rapidly form chains which

111

 

then aggregate into thick columns [Ha1590, Chen92, Mart92, Gind92]. Magnetic
ferrofluids form similar structures of chains and thick columns [Wang94]. The de-
pendence of the periodicity of columns on the sample thickness has been recently
studied using magnetite colloids [Wang94]. Also many numerical simulations us-
ing molecular dynamics (MD) [Hass89, Klin89, Kusa90, Bonn92, TaoR94] have

been carried out to find the orientational ordering in ER liquid crystals.

Even though a few MD simulation results have been published for ferroelectric
dipole liquids [WeiD92a, WeiD92b], almost no MD simulation results are available
for the magnetic ferrofluid systems except a few simulations using the Monte Carlo
technique[WeisQ3, Stev94]. Although ER fluids and magnetic ferrofluids have
many common properties there is one fundamental difference. In ER fluids, the
dipoles are induced by applied electric field and they are always aligned with the
field. Consequently, the dipole moments vanish in the absence of a local electric
field and therefore do not show any interesting behaviour in that situation. In
contrast to this, the magnetic particles in ferrofluids have their own pre-assigned
magnetic moments. This means that in general case they are not necessarily
aligned with the external field, and we have to consider the rotations of these
particles as well as their positions. Also, we should expect that these fluids have

some non-trivial structures even in the absence of the external field.

In this Chapter, I will discuss the early stages of colloidal aggregation in a
ferrofluid and the stability of the ring structure of these particles. In an ER
fluid it is very unlikely to observe such a ring structure for the reasons I give in
the previous paragraph. As I will show in later sections, a ring is the most stable

configuration in zero-field, where it is established by the dipole-dipole interactions.

In an ER system such interactions only occur upon applying an external field.
112

 

 

 

 

 

 

 

 

 

 

 

 

In Sec. [5.2] I will describe the theoretical model used to simulate the mag-
netic ferrofluid. The interaction potentials for these colloidal particles will be
introduced. Since the magnetic colloidal particles have permanent magnetic mo-
ments pointing in a particular direction, we need to treat their rotational degrees
of freedom properly, as well as translational. Sec. [5.3] will be devoted to the
techniques involved in doing MD simulations for rigid bodies, especially spheri-
cal tops. Since the standard equations of motions for the Euler angles which are
used to represent the rotation of rigid bodies, have a well-known singularity, the
quaternion formalism is introduced in Sec. [5.4] as a better alternative. Before we
use all these techniques and start the MD simulations, in Sec. [5.5] the stability of
the chains and rings made of magnetic dipoles is discussed under various physical
conditions with static or energy minimization method. The dynamical response
of the magnetic dipole structures is the subject of Sec. [5.6]. The phase diagram
of the ring configuration will also be presented. The ongoing projects in ferrofluid
will be briefly described in Sec. [5.7] and in Sec. [5.8] I summarize this chapter

and draw conclusions.

This chapter contains material of the following two separately published

manuscripts:

0 [Jund94] P. Jund, S.G. Kim, D. Tomanek and J. Hetherington, Stability of
Magnetic Dipole Structures in Ferrofluids, preprint (1994).

0 [Kim894c] S.G. Kim, P. Jund, D. Tomanek and J. Hetherington, Structure

Formation in Ferrofluids, in preparation (1994).

113

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

5.2 Modeling of the ferrofluid

The physical system we attempt to simulate in this study is very much related
to the ferrofluid system used in Ref. [Wahg94]. The sample used by Wang et
al is a magnetic colloid made of magnetite particles coated with oleic acid and
suspended in n-eicosane. The volume fraction of magnetic particles is 12% with
mean diameter of particles 89 A. The thickness of the coated surfactant layer is

about 20 A.

Based on these experimental data, we choose 100 A for the diameter of parti-
cles. From the density of the magnetite, 5.2 g/ cc, and their size we also set the
mass of our magnetic dipoles to be 1.64 x 106 amu. Using the magnetization of
magnetite and the volume of the particles, we can estimate the magnetic moment
of these particles by assuming that they have a single magnetic domain. We take
2.1 X 104mg, with pa being Bohr magneton, to be the amplitude of the magnetic

moments of all particles.

In this calculation, the ferrofluid particles are modeled by soft spheres of di-
ameter a and carrying a magnetic dipole moment ii = pail. We take no to be
constant for all dipoles and it is the unit vector along the magnetic dipole moment.

The potential energy of the system of these dipoles can be written as

U=2 [u.(i)+2:u(m] on

be

where the pair potential between particles i and j is

u(tj) = 1109(5).) + usc(‘ij). (5.2)

114

 

 

The first term in Eq. (5.1) is the field energy
u3(i) = —;z.- . r3 (5-3)

gained by each particle when an external magnetic field, E is applied.

The magnetic interaction between the particles is

upo(ij) = it.-fl.-/r?,-—3(fl.~fi-.-)(fi.--r':-.-)/r?.-

= 573 [iii 111' — 301; 'fiin‘j 4.3)] (5.4)

51'
which is the usual dipole-dipole interaction[Jack75]. Here ii,- and F.- are the dipole
moment and position vector of particle i, respectively. Obviously, r}, = 1'",- - 1"},
and Ti,’ is the magnitude of 1"},- and 5.3 is the unit vector along 1%,. Using the
magnetic moments and diameter of the particles the energy scale of the dipole-

dipole interaction is

2
”5% g 2.368 x 10‘2eV. (5.5)

There has been a significant amount of efforts to measure the pair interaction
potential between the colloidal particles [Cald94, Croc94]. In this study we use

following general form for soft-core potential

u$605) = 69(1's/0') (5.6)

with

 

 

<I>(r) = a [exp (J. p.10) - exp (‘1. go” ’ P1 < P2- (5-7)

115

 

 

This functional form is similar to the one suggested by Tejero et al. [Teje94]
except the 1/r factor which is in their potential. Here 6 fixes the energy scale of

this potential and o is the distance where
1150(0) = 0 (5.8)

and defines the length scale of the soft-core. The three parameters (a,p1,p2) of
Eq. (5.7) are seen to determine the steepness of the repulsion (p1), the range of

the attraction (p2), and the well depth (a). If we fix the well depth to be s,
@(m) = -—1 (5'9)
with re denoting the position of the minimum of <I>(r), viz.
<I>’(ro) = o. (5.10)

This condition will fix a = a(p1, p2) for given p2 and pg

1
(&)m/(pz-m) .. (fl)m/(m-m) °

PI PI

 

a:

(5.11)

To simulate the effect of surfactant layer in ferrofluid, we choose the depth of the
potential well to be about 10% of the dipole-dipole interaction, or 2.0 x 10‘3 eV
and the range of two exponential terms to be p1 = 2.5 A, p2 = 5.0 A. These
parameters gives a = 4 and generates the potential which decays to 10% of its
depth. at the distance of surfactant layer thickness, 20 A, away from the particle
diameter o, 100 A. Fig. 5.1 is a plot of uso(r) with this set of parameters. The
maximum of the plot is 2.0 x 10"2 eV which is the typical value for the dipole-

dipole interaction.

116

 

 

200

160

Energy [10“ 0V]
an
o

O
O

 

 

 

 

 

ITI I I r rl—TI I1 l—rlrlTj {Tl [FT [j
- -
I- -I
I- d
1— —(.
-_ -I
- -
I- ..
- -
— —.(,
- -
I- -
b u
I- -I
I— —.
I- c-
n ..
V «II
. \ -
I. W -I
I- .
h '1
- I J L. I | -

I I I LI I I I_I I 7 I_I L LI I_I ,J I I I

 

 

 

100

110 120 130 140

Distance [A]

Figure 5.1: Soft-core potential for a ferrofluid ' ° '
particle, use r , wrth e = 2.0 1 3
eV,p1=2.5Aandp3=5.0A. U x 0

117

 

 

 

 

 

 

For completeness, I summarize the units and conversion factors between these
units in Table 5.1. All the units with tilde are the units used in our numerical
calculations. In Table 5.2, the actual physical quantities used and the typical

values we need to deal with in this study are listed.

The dipole-dipole interaction in Eq. (5.4) has two terms competing each other.
The first term, it,- 42,-, tends to make the dipoles anti-parallel each other and favors
anti-ferromagnetic arrangement. On the contrary, the second term, —3(fi.- -f,-J-)(fij .
1%,), tends to make the dipoles align themselves and favors a head-to-tail “chain”
structure. Since the second term has a three times bigger prefactor, it is easy
to understand why the chain structure is most commonly observed in ferrofluids
and other strongly interacting dipole fluids. Suppose we have two parallel dipoles
whose positions are fixed along z-axis and make angle 6 with z-axis as in Fig. 5.2.
The first term in Eq. (5.4) is fixed to unity while the second term changes as we
vary the angle 9. When 9 is small the dipole-dipole interaction is attractive, but

above a certain angle, 9c, the interaction becomes repulsive. This angle satisfies
1 — 3cos2 9,, = 0 (5°12)

which yields 9., = cos‘1 1/3 or 547°. This effect will be seen again when we

consider the breaking of the rings of these particles.

The formulation is exactly the same as for electrorheological fluids [TaoR94].
There the particles have electric dipole moments, 5,, and are subject to an external
electric field, E. One simply needs to replace all [[38 with 131’s and E with E. But
in this latter case, the electric dipole moments are not fixed, but rather induced

by the applied field.

118

 

 

 

Table 5.1: Units and conversion factors of physical quantities related to the present
model of ferrofluids.

 

 

 

 

 

 

 

 

 

 

 

 

 

Physical quantity Symbol Value
Length A 100 A
Mass amu 106 amu
Energy ell 10" eV
Time = AW 1.02 x 10‘7 S
Magnetic moment a}; 1.3647 x 103 p3
Magnetic Field é 12.6592 G
Temperature K 1.16049 K
Dipole-dipole energy [32 /A3 ev
Dipole-field energy I53 C eV
Thermal energy kBK ell

 

 

 

 

 

 

119

 

 

 

Table 5.2: Physical parameters and typical values of some quantities used in
ferrofluid study.

 

 

 

 

 

 

 

 

 

 

 

 

 

Physical quantity Symbol Value
Length 0 1.0 A
Mass m 1.64 amu
Magnetic Field B 10.0 G
Magnetic moment no 15.388 p'g
Dipole-dipole energy pig/0'3 236.8 eV
Dipole-field energy [‘03 153.88 ell
Soft-core energy 6 20.0 eV
Soft-core repulsion range p1 2.5 A
Soft-core attractive range p2 5.0 A

 

 

 

 

 

 

 

120

 

 

 

 

 

 

Figure 5.2: Two dipoles on :s-axis.

121

 

5.3 Molecular dynamics of spherical tops

Since we are modeling ferrofluid particles by spheres carrying a dipole moment
each particle has not only positional coordinates, but also orientational coordi-
nates. The former set of coordinates is very common in any molecular dynamics
simulation, but the latter requires some careful considerations[Alle90]. The orien-
tation of a rigid body specifies the relation between an axis system fixed in space
and one fixed in that body. Usually the ‘principal’ body-fixed system in which the
inertia tensor is diagonal is chosen. In our case we choose the third axis along the
dipole moment and the other two axes on the sphere’s equator plane. Any unit
vector 5' may be expressed in terms of its components in body-fixed or space-fixed
frames: we thus use the notation é” and 5', respectively. The rotation matrix A

relates these components by
at = A . 51'. (5-13)

The nine components of the rotation matrix are the direction cosines of the body-
fixed axis vectors in the space-fixed frame, and they completely define the molecu-
lar orientation. In fact there is substantial redundancy in this formula: only three
independent quantities (generalized coordinates) are needed to define A. These
are generally taken to be the Euler angles (MW in a suitable convention (See Fig.

5.3 and [Gold80]). Then,

A = . .
- ' ' ' ¢+cos¢cos08in¢ smflsmip 5.14
f::f;::n¢¢ -s;riin¢¢c::89;::;/J 1), .8 131:“; :isn 1p + cos ¢ cos 9 cos it sin 0 cos ¢ ( )
sin 45 sin 0 — cos 4’ Sin 9 cos 9

122

 

 

Figure 5.3: Definition of Euler angles.

123

 

 

 

 

 

Clearly, if E is a vector fixed in the particle body-frame (such as the dipole moment
vector), then é” will not change with time. In spacefixed coordinates, though, the
components of 3' will vary. This is a special case of the general equation linking

time derivatives in the two coordinate systems[Gold80]
é=é+an=aUdfi aw)

The time derivative of the angular velocity vector I3 is dictated by the torque
ii acting on the dipole. Although the torque is most easily evaluated in space-
fixed axes, it is most convenient to make the connection with 6, via the angular

momentum l: We use the Newtonian equation
I", ___ 1.1:. (5.16)

and obtain, after applying Eq. (5.15), the expression

i;=ii’+t3' Xi; (5'17)
01'
a+axr=s. mm)

Eq. (5.18) is thus the appropriate form of the Newtonian equation of motion

relative to the body axes. The ith component of Eq. (5.18) can be written as
ii + Egjgwjlh = 11;, (5.19)

(Repeated indices imply summation.)

where the “body” superscript has been dropped. If now the body axes are taken

124

 

 

 

 

 

as the principal axes relative to the reference point, then the angular momentum

components are l.- = Igwg. Equation (5.19) then takes the form
1.1.ng + 6,1"):ijka = 11;, (5.20)

since the principal moments of inertia are, of course, time independent. For
spherical tops, such as our spherical dipole moments, all three principal moments

of inertia have same values I and the cross product term in Eq. (5.20) yields
0‘6 = ”/1. (5.21)

Conversion between the space-fixed and the body-fixed system is handled by the
analogues of Eq. (5.13), i.e.

a“ = A .11" (5.22)
a" = A“ .sbzAT-o'b (5.23)
since the inverse of A is its transpose. To complete the picture, we need an

equation of motion for the particle orientation itself, i.e. for its rotation matrix

A. We may write down the equations of motion of the Euler angles [Gold80]

 

. _ ‘sin¢cos9 ,COB¢C°39 I

4’ — —w: sin9 " sin9 ’

é = w;cos¢+w;sin¢

. _ ,sin¢_ ,cosdi

I“ — w’sin9 w”sin9° (5.24)

These equations would apply to each dipole separately to obtain its own new
rotation matrix A. In principle, Eqs. (5.14) with (5.21)—— (5.24), may be solved

in a step-by-step fashion, just as we deal with translational equations of motion.

They, however, suffer from a serious drawback. The presence of the sin 9 terms

125

 

 

 

ll

 

 

 

 

in Eqs. (5.24) means that a divergence occurs whenever 9 approaches 0 or 11'.
The molecular motion is of course, unaffected when this occurs but, because of
our choice of axes, the angles ()5 and 1]) become degenerate. Thus the equations
of motion are unsatisfactory when written in this form. One way to cope with
this would be to reduce the time step, so as to deal with the more rapidly varying
quantities, whenever any molecule approached the critical values 9 z 0 or 1r.
This would be very expensive and awkward. A slightly more satisfactory solution
was employed by Barojas, Levesque, and Quentrec[Baro73]: two alternative sets
of space-fixed axes for each molecule were used, and whenever the angle in one

system approached 0 or 1r, a switch over to the other set was made.

In recent years, a much more elegant and straightforward solution to the prob-
lem of divergence in the orientational equations of motion has been proposed.
Recognizing that singularity-free equations could not be obtained in terms of three
independent variables, Evans [Evan77a] suggested the use of four quaternion pa-
rameters as generalized coordinates. Quaternions fulfil the requirements of having
well-behaved equations of motion. The four quaternions are linked by one alge-
braic equation, so there is just one ‘redundant’ variable. The basic simulation
algorithm has been described by Evans and Murad[Evan77b]. In the following
section, I will briefly summarize the theory of quaternion parameters adapted in

our molecular dynamics simulations.

5.4 Quaternion molecular dynamics

The procedures to utilize the quaternion parameters in molecular dynamics has

been summarized beautifully by M.P. Allen [Alle84].

126

 

 

 

 

 

 

A quaternion Q is a set of four scalar quantities, or a combination of a scalar

and a three-component vector

Q = (QOi91i92aq3) = (90,4). (525)

The “quaternion product” of two quaternions P and Q is defined via the combi-

nation of the usual scalar product 13' - (j' and vector product 15' x q" as

P*Q=(poQo—i-§3poci'+qoi+i><<f)- (5-26)

One thing to note is that there is a sign mistake in Allen’s definition [ Eq. (B2) in

[Alle84]] which could confuse many readers. As I verified in Eq. [Evan77a] (5.26)
is correct. Of course, the scalar product of two quaternions is still defined in the

usual way
P-Q=poQo+i-é'- (5-27)
The quaternion product is associative
P*Q*R=P*(Q*R)=(P*Q)*R, (5.28)
but not commutative

P*Q aé Q*P. (5.29)

The conjugate of Q is defined as

Q = (90: “-6): (5'30)

so that

6* Q = as Q = (IQI’fi), (5.31)

 

 

 

 

 

 

 

where IQI2 = Q - Q is the norm of Q. The quaternions of interest here satisfy the

constraint
IQI2 = <13 + 9i + 43 + 43 = 1- (532)

Therefore Q is the same as Q“, the inverse of Q, since the quaternion (1, 0) acts
as an identity element in the quaternion product. The way in which a quaternion
with unit norm may represent the orientation of a rigid body is discussed in Sec.
4-5 of [Gold80]. For example, the transformation between body-fixed vectors 1:6

and space-fixed vectors 1‘” may be expressed as
R” = (sumo, (5.33)

where R” = (0,1‘”) and R‘ = (0,13).
The equations of rotational motion for each particle are advantageously ex-
pressed in quaternion form as

.1.

20' s: Q, (5.34)

. 1
Q= §Q*nb=

where 0" = (0,6”) and 0‘ = (0,132“). Here, :3” is the body-fixed, «3' the space-
fixed, angular velocity. For a spherical top, following Eq. (5.21), the angular

velocities evolve according to
(I5 = Nb/I or fl‘ = N‘/I, (5.35)
where I is the moment of inertia, N” = (0,175), N‘ = (0,1'1"), and ii" and it" are

respectively the torque in the body-fixed and space-fixed frame.

These first-order differential equations can be integrated in time by various

standard methods[Alle90]. But calculating the torque directly from the molecular
128

 

 

orientations is not simple when many dipole moments are present. It would be
nice if we could express the torques in terms of a “quaternion gradient” of the
potential energy. At this point, we will eliminate the angular velocities from these

equations. Differentiating Eq. (5.34) with respect to time we get

- 1 . b 1 . b

Q=§Q*n +§Q*n. (5.36)
Using Eq. (5.34) we note

1 b " ‘

5n = Q i. Q. (5.37)

Taking derivatives of both side of Eq. (5.31) and applying the quaternion product

of Q* from the left, we get

Q .. Q * Q = —IQI2Q- (5-38)
Combining these results we arrive at the second order differential equation

” ‘ 2 1 5 39

Q=-IQ|Q+ZN/I, (- )
where

N = 2Q *N" = 2N‘ * Q. (5-40)

The elegant feature of this procedure is that N can be expressed directly in terms

of quaternion parameters[Alle84]:
N = —(¢9U/6qo,6U/6q1,6U/0q2,6U/6q3). (5.41)

Here U is the potential energy of the system as a function of the quaternion

parameters of all the particles, and we differentiate with respect to the quater-

nions of the particle under consideration. These derivatives should really be taken

129

 

subject to the constraint, Eq. (5.32), but then it would be impractical to use Eq.
(5.41) to solve the equation of motion. Fortunately,-however, it is possible to treat
the quaternions as independent parameters, and take unconstrained derivatives,
if the consequent motion off the constant-IQI2 hypersurface is projected out. This

correction should be performed at every step in the following procedure [Alle84]:

[Step 1] After calculating N from Eq. (5.41), we subtract the quantity
(Q - N)Q, which yields

N' = N — (Q - N)Q. (5.42)
[Step 2] Using N’ obtained in [Step 1], calculate N” or N‘ from Eq. (5.40)

Nb: Q*N’ or N'=%N'*Q. (5.43)

NIH

These quantities are used in Eq. (5.35) to get new II" or fr.

[Step 3] Advance (2" or 0' one time step to get their new values, using

the result of [Step 2].

[Step 4] After calculating Q from Eq. (5.34), the quantity (Q - Q)Q is

subtracted, yielding r
Q’=Q-(Q-Q)Q- <5-44)

[Step 5] Advance Q one time step using Q’ obtained in [Step 4] to get

new quaternion parameters.

[Step 6] Renormalize this new Q to have unit norm.

130

 

 

Note that after these corrections N’ and Q’ are perpendicular to Q, yielding

Q-N’ _ o

Q-Q' = 0. (5.45)
so that

Q - Q’ = -IQI’. (5.45)

In the above procedure, the [Step 1], [Step 4] and [Step 6] are the critical steps
to keep our quaternion parameters on the constant-IQI2 hypersurface. We can
understand this more clearly by expanding Q over the time domain. Suppose

that we have normalized Q at time t,

|Q(t)l’ = 1. (5.47)

At time t + 6t,
Q(t + 5) = Q(t) + Q’(t)5t + $61059 + 0(6t3). (5.48)

The the norm of Q becomes

IQ(t + 6t)!” = IQ(t)I2 + 2Q(t) - Q’(t)at + IQ’(t)I’6t’ + Q - Q’stz + 0(6t3)
= 1 + 0(5t3). (5.49)
Thus previous correction procedures guarantee us conservation of the norm of Q
up to the second order in it without renormalizing.

Another important reason to implement this correction comes from the self-

consistency of the formalism. By definition, the 0-th (or “scalar”) component [See

131

 

 

 

 

definition in Eq. (5.25)] of N” or N‘ should vanish. From Eqs. (5.42) and (5.43)

N' = %IN*Q-(Q-N)Q*Q]

= éINaeQ—(Q-NXLO'N. (5-50)

Recalling the definition of Q, we can easily show that the 0-th component of
N‘ vanishes. Futhermore we need to note that the correction [Step 1] not only
eliminates the O-th component of N' but also it is the only modification it makes
as far as N‘ is concerned. Thus we could skip the correction [Step 1] and do the
[Step 2] and simply discard the 0-th component of N” or N‘. A similar argument
can be used to show that the correction [Step 4] eliminates the 0-th component

of (1" or 0‘ if we have to calculate them by inverting Eq. (5.34) to
0" = 2Q * Q’ or 0' = 2Q'*Q. (5.51)

In the Euler angle convention of Fig. 5.3 and Eq. (5.14), it is most convenient

to define
qo = cos $9 cos $035 + 1/1)
q1 = sin %9 cos $00 - 111)
q; = sin %9 sin %(¢ - 1/1)
q; = cos -;-9 sin $0” + 1,0). (5.52)

Then the rotation matrix becomes

2(9192 - 9093) 95 "' 9i ‘l’ 93 " 95 2(9293 + 9091)

I 93 + 9? - 9% - 9% 2(9192 + 9093) 2(9193 - 9692)
A = , (5.53)
2(9193 + 9092) 2(9293 - 9091) 93 - 9i - 9% + 93

which is obtained by converting Eq. (5.33) into matrix equation. Similarly the

132

 

 

 

 

 

 

equation of motion for quaternions for each particle, Eq. (5.34), can also be

converted to the simple matrix equation

90 90 -91 -92 —9a 0

91 = _1_ 91 90 _q3 q2 w: (5 54)
92 2 92 93 90 “91 w; '

93 93 —92 91 90 01,

As we expected, the equations of motion, Eq. (5.54) with Eq. (5.21), contains no
unpleasant singularities, if the transformation matrix (5.53) is used, to transform

between space-fixed and body-fixed coordinates.

5.5 Rings and chains

Fig. 5.4 depicts the structure of a chain and ring of 10 magnetic dipoles. The
directions of magnetic dipoles are represented by two hemispheres. Blue (dark)
and orange (light) color represents the north and south pole of the dipoles, re-

spectively.

Due to the dipole-dipole interaction, described in Eq. (5.4), the magnetic
dipoles want to align. For a small number of dipoles in the absence of exter-
nal field, a chain structure is the most stable configuration. As the number of
dipoles increases, ring configurations become more favorable over chains, since
the increase of the energy due to bending of the aligned dipoles will be more
than compensated by joining two ends. The potential energy per particle for the
chain (open diamond) and ring (solid square) configuration as a function of the
number of the magnetic dipoles is plotted in Fig. 5.5. Each energy is obtained

as a result of complete structure optimization using the conjugate gradient (CG)

method [Pres86]. There is no external field applied. We note that for the number
133

 

 

 

 

 

 

Figure 5.4: 10 ferromagnetic spheres in a (a) chain and (b) ring structure. The
directions of magnetic dipoles are represented by two hemispheres. Blue (dark)

and orange (light) color represents the north and south pole of the dipoles, re-
spectively.

134

 

 

 

 

 

 

IIIIIIIIIIIIIIIIIIITIII

-300 I '-

E(N)/N [10" eV]

-600

 

 

 

_600IIIIIIIIIIIIIJLIJIIILIII

0 5 10 15 20 25
Number of particles N

Figure 5.5: Energy per particle of a chain (open diamond) and a ring (solid square)
0f magnetic dipoles as a function of the number of particles in the absence of
external field. The solid and dotted lines are the predictions from Eq. (5.61) and
139- (5.62), respectively, with a,=1.20 and ac=1.18.

 

 

 

 

 

 

of dipoles, N is bigger than 4, rings are always more stable than chains. The more
dipoles are in the ring, the smaller is the bending angle and energy. It also should
be noted that the difference in energy is maximum around N 2 9 and is getting
smaller as N increases. This behaviour can be understood easily with very simple
arguments. Consider a N dipoles ring as in Fig. 5.6. Since we have N dipoles,

the angle between two neighboring dipoles is
9 = 21r/N (5.55)

From the figure we can easily see that neighboring dipole 1 and 2
it: #12 = cos 9

[11 ° F12 = C08 59

= [22 ° P12 (5.56)

and the dipole-dipole interaction between them is

2 3 9
U12(9)=-2% (312;) (5.57)

We know that the distance between the dipoles d will change for each different N
to minimize the energy. In our simple argument we assume that it is independent
of particle number N and set to their diameter 0. Furthermore we will neglect the
details of the higher-order neighbor interactions and simply say that the energy

for N particle system is
Um(N) z a,NU12(9) (5.58)

where a, is the constants to be determined and independent of N and contains

the information on the soft-core potential. Thus the energy per particle becomes

Ufin5(N) ~ pg 3 + cos(21r/N)
T _ —2a,— 4 ] (5.59)

136

 

 

 

 

 

Figure 5.6: A vector representation of ring of 10 magnetic dipoles. The directions
of magnetic dipoles are represented by the arrows.

 

 

 

 

 

 

 

and with similar arguments for linear chain

Uchain(N) pgN—I
N _ —2aca3 N
_ _ s( -1)
._ 211303 1 N (5.60)

By matching the values from Eq. (5.59) and Eq. (5.60) to the ones obtained from
full energy-optimized calculation at large N, say 25, we can set the values of a, and
ac. The solid and dotted lines in Fig. 5.5 are the plots of Eq. (5.59) and Eq. (5.60)
with a,=1.20 and ac=1.18. The analytic curves represent the numerical data very
well except a little deviation in small N region. This is a remarkable agreement,
considering how simple and crude our arguments are. We may also notice that it
even predicts correct N value to make the ring more stable, which is 4. Also from

these analytic expression, we can get the scaling behaviour for large N

«21

, 2

_ . . . . 5.61
N + ) < )
Thus the energy per particle for the ring decays as 1/ N 2 while the energy per

particle for the chain decays as 1/N.

.9 Now let us apply an external magnetic field to these systems. The chains will
obviously align with the external field to minimize the energy. The behaviour of
the ring is not very trivial when subject to an external field. They would like
to align with the field but then half of them have to be “anti-parallel” to the
field. This increases the total energy , resulting no net energy gain if they remain
a perfect circle. The rings can gain energy slightly by distorting their structure
from the perfect circle. Fig. 5.7(a) is a 20 particle ring forced to remain parallel

to the field. The external field is applied horizontally and is represented by an

arrow in the figure. We used the same coloring convention for the dipoles and the

138

 

 

 

 

 

 

field. The dipoles in Fig. 5.7(a) are allowed to adjust their positions in the plane
of the ring. The ring makes an “egg-shape” which has a flatter side and a pointed
tip. Having the flatter side makes the ring accommodate “more” dipoles parallel
to the field. The pointed tip, however, costs a lot of energy because of the sharp
turn of the dipole direction and becomes very vulnerable as the field increases.
Above a certain critical magnetic field, this tip becomes too sharp and the ring
can not support the stress any more. The ring breaks at this tip and becomes a

chain or fragments to little pieces, depending on the surrounding medium.

It turned out, however, that the rings can stand much higher fields without
breaking by being perpendicular to the field, as shown in Fig. 5.7(b). In the
figure the field is pointing out of the paper, as represented by the head of the
cylindrical arrow. It should be noted that their positions remain on a perfect
circle, yet the dipoles show a tendency to align with the field. In figure, you
see more blue (dark) colors for each sphere meaning that they are facing toward
you. The decreasing dipole-dipole interaction within the ring is compensated by
the field-dipole interactions. Above a certain critical field, these dipoles rotate
sufficiently out of the ring plane to become nearly parallel to each other, and their
mutual interaction vanishes. At this moment the ring flies apart into individual

particles.

Therefore in the following discussions I will focus only on the chains parallel
to the external field and rings perpendicular to the external field. The energies
per particle of a chain and a ring of 10 dipoles are plotted in Fig. 5.8. As seen
in Fig. 5.5, at zero field, the ring structure has the lower energy. As the field

increases, the energy per dipole of the chain decreases linearly in the external

field. In contrast to this, the energy per dipole of the ring structure remains almost
139

 

 

 

 

 

Figure 5.7: Rings of 20 magnetic dipoles in an external field which is (a) parallel
(300 G) and (b) perpendicular (500 G) to the ring plane. The external field is
represented by the arrows, and has the same coloring convention as the magnetic
dipoles, used in Fig. 5.4.

140

 

 

 

 

 

 

 

 

 

 

 

IIIIITIIIIITIIIIIITIIlllq

 

 

 

1
_2000 MIIIIJIIJ’ILLIIIJJJIJJII

200 400 600 800 1000
Magnetic field B [Gauss]

 

Figure 5.8: Energy per particle of a chain (open diamonds) and a ring (solid
squares) of 10 magnetic dipoles as a function of the external field.

141

 

 

 

constant at low fields and decreases very slowly at higher fields. At a relatively
- low field value, about 37 G for N = 10, these two curves cross and the chain
structure becomes energetically more favorable. Once the ring structure is formed,
however, the ring does not break into a chain at this external field. The ring is
in a local energy minimum in the multi-dimensional coordinates space and there
is an activation barrier to overcome. We can go over this barrier dynamically by
pumping energy into the system and increase the kinetic energies of the particles.
When the added total energy added exceeds the height of the barrier, the system
will find the saddle point in phase space and move into the lower minimum. This
is the subject I discuss in following sections using the molecular dynamics study.
The alternative way of finding the lower (local) energy minimum is to lower the
barrier by changing the external field. When we increase the external field, we are
changing the contour of the energy surface in multi-dimensional coordinates space.
Increasing the magnetic field will lower the energy barrier for a certain geometry,
which will result in a saddle point pathway leading from the higher local energy
minimum of ring structure to the lower energy minimum of chain structure. For
10 dipoles ring, this energy barrier is removed at external fields exceeding 750 G,

as shown in Fig. 5.8.

Therefore, from Fig. 5.8 and the preceeding arguments we can define two
critical fields for the rings of a given number of dipoles. The lower critical field,
301 , is the field above which a chain becomes energetically more favorable. The
upper critical field, 3.32, is the field at which the energy barrier between the chain
and ring structure disappears. Fig. 5.9 is the plot of Bel for a given magnetic field.

The results obtained from conjugate gradient method (squares) are compared with

142

 

 

 

 

the analytic predictions (curve) using Eq. (5.59) and Eq. (5.60)

U____chain( N) _U—_ring( N)

 

 

floBcl = N
2 2"—° cap —a.) +N — °+a.c°s(2":N) - 1]. (5.62)
And for large N
a¢_ am”
#03012 2— a:[(—a, ac) +— N— 2N2 + ~ - ] . (5.63)

This analytic expression predicts correct 1/N scaling behaviour for large N and
even predicts correct N to have maximum Ben which is 9. Fig. 5.10 is the plot of
BC: for a given magnetic field obtained from conjugate-gradient energy optimiza-
tion. BC: is related to the height E‘ of the barrier between the minimum energy

configuration for the ring and that of chain. We can estimate E" using
E‘ N FOBc‘b (5.64)

as the amount of energy needed to break a ring into a chain or another lower
symmetry structure in the absence of the external field. One thing to note is that
these estimates are the result of straight energy minimization calculations with
absolutely no dynamics. In next Section we will consider the stability of dipole

rings with the full dynamics.
5.6 Molecular dynamics study

We studied the stability of the magnetic dipoles by considering the ground state
energies of ring and chain structure. This study, however, only gives the estimation

of the energy barrier in zero field. If we want to know the stability of the ring for

143

 

 

 

 

50llTTIFlFIIIIIIIIII Ill

30

IIlIIIIlIIII

  
  

 

 

20

Magnetic field [Gauss]

10

[Llllllll

 

 

IIIIITTIFTIIIIIIIIIIIII

 

IIIJILIIILLIIIIIIIIJJ

0 5 10 15 20
Number of particles N

 

N
0'

Figure 5.9: Lower (Be-1) critical field for magnetic dipole rings as a function of the
number of particles. Results from energy minimization (squares) are compared to
analytic expression (curve) from Eq. (5.62).

144

 

 

 

1000 UFIMIMIII'IITI‘IT'IWITIMTTI

 

800

600

Magnetic field [Gauss]

200

TsllrlsrlsrllllITllitrl

IIIIIIIIIIIIIIIIIIII

 

 

 

IIILIILIIIIIIJIIIIII

5 1O 15 20 25
Number of particles N

I

0

Figure 5.10: Upper critical field (Bc2) for magnetic dipole rings as a function of
the number of particles.

145

 

 

 

 

 

a given field, we need to do the dynamics study. For MD simulation we integrate
the equations of motion, such as Eq. 4.29, with fourth-order Runge-Kutta method
[Pre586]. The equations of motion include the quaternion parameters to represent
the rotation of the dipoles. We use an adaptive step-size control in which the time
step 9t varies according to the moves of particles. Two criteria, 6rd and 5702 are
specified. If the largest position change among all particles, firm“, is smaller than
6rd, the time step «it in the subsequent integration will be increased. When 57:5“,
is bigger than 5mg, the time step 5t decreases. Otherwise we keep the current time
step 5t. All of the MD results in this chapter are generated with the initial time

step 6t = 5.0 x 10‘11 S, 6rd = 0.8 A and 5rd = 1.2 A.

For a given field, we first obtain the optimized configuration using the
conjugate-gradient method. We then give a small kinetic energy corresponding to
a temperature of 1 K. The velocities and angular velocities are generated randomly
according to the Maxwell-Boltzmann velocity distribution of given temperature.
We use the step-wise microcanonical MD simulation meaning that we pump en-
ergy into the system at one time and keep the system isolated for a sufficiently
long time. The system will explore the multi-dimensional phase space available
within that given total energy. If the barrier is lower than the total energy the
system will eventually find the pathway leading to other local minimum. After
a suficiently long time, if ring is still intact, it usually means that the barrier is
higher than the total amount of energy pumped in. After then we increase the
total energy of the system by heating up the particles (scaling of velocities) and
repeat above simulation steps until we have a broken ring. The amount of energy
we pump in at each step will determine the resolution of energy scale of our study.

In Fig. 5.11, we plot the variation of the energy barrier, E‘, as a function of the

146

 

 

 

 

 

applied field, B. To generate the data we run MD simulations for an isolated
system (microcanonically) at each given total energy for 1.0 X 10'6 S. We add
3.0 x 10"3 eV of energy into kinetic energy if the ring is not broken until that

time.

E‘ is the amount of total energy you need to put in to break a ring of ferrofluid
particles at given magnetic field. In the region below this curve, once you have
a ring of ferrofluid particles, you will have a stable ring. When the total energy
difference from the optimized configuration is more than this value, the ring will
disintegrate into chains or smaller pieces. According to the “equipartition theo-
rem” [Reif65a], the average kinetic energy, I? of a ferrofluid system which has 6

degrees of freedom per particle (3 translational and 3 rotational) is
- 1
K = 6 - 2kBT = 3k3T. (5.65)

Fig. 5.12 is the plot of the estimation of melting temperature TX, as a function of
magnetic field. Ti! is calculated from Eq. (5.65) using the average kinetic energy

of the system when the ring breaks.

5.7 Ongoing work

The ferrofluid as well as an electrorheological fluid is a system of dipole parti-
cles immersed in a viscous fluid. It is well known that these sufficiently small
macroscopic particles exhibit a random type of motion, so—called “Brownian mo-
tion” [Reif65b]. If we consider the effect of this Brownian motion, the governing

equation of motion becomes the “Langevin equation”

mid; = y: _ 10 + F’(t). (5.00)

 

 

 

 

 

BOOIIrIIIrIIIIIfIIFII

 

 

500 __ --E

E i

400'-' -*

if : :
T Z I
c 300— -'
1'2. I 3
a I 3
200:7' -:
1005- —E

o ' I I I 1 l 1 4 LI J L 1 1 1 I :

 

 

O

200 400 600 800
Magnetic field B [Gauss]

Figure 5.11: Energy barrier E" of ferrofluid rings as a function of applied field B.

148

 

 

 

ZOOIrIrIfiIIIIIIfil'ITIF

.

q

q

d

-

d

a" -

I_I _,

.3- t

u

q

-

—|

q

.

q

| | '
01L11_4LmLIILL

 

 

 

0 200 400 600 800
Magnetic field B [Gauss]

Figure 5.12i Estimation of melting temperature T5, of ferrofluid rings as a function
of applied field B.

149

 

.7: is the force that we have without the Brownian motion, such as the one from
the dipole-dipole interaction. The frictional coefficient 7 is given by “Stoke’s law”
[Page52, Joos58]

7 = 31rna' (5.67)

where 17 is the viscosity of the surrounding fluid and a is the diameter of the
particles. The magnitude of the Brownian force depends on the temperature T
and the viscosity 1) of the fluid. To simulate the effect of Brownian motion we
generate the random force F’ (t) which has the Gaussian distribution with zero-

mean and the width

 

fl = \/61rlc3To'n/r (5.68)

where r is an adjustable parameter and has the dimension of a time [TaoR94].
With given values of the state parameters, we can tune r so that it generates
ensemble of systems with the average temperature matching the heat bath or fluid
temperature. We have very preliminary MD simulation results with this model of
Brownian motion and it will be part of our next paper. Also, we are simulating
the aggregation of many of these ferrofluid particles into chains and branched
structures. Fig. 5.13 shows one of our preliminary results on the aggregation of
64 magnetic dipoles in the absence of the field. At the start of the simulation,
64 magnetic dipoles are perturbed randomly by 10% of the lattice constant, 173
A, from their 4 x 4 x 4 cubic lattice structure and only the top layer is seen
clearly. After about 1.0 x 10"8 S, the particles aggregate into two rings and one
giant “snake”. One exciting fact is that this aggregation was done in zero-field.

Many simulations and experiments have been performed in rather high-field and

150

 

 

 

Figure 5.13: Aggregation of 64 magnetic dipoles. (a) Initial configuration and (b)
final configuration. In (a) 64 magenetic dipoles are perturbed randomly by 10%
of the lattice constant, 173 A, from their 4 x 4 x 4 cubic lattice structure and only
the top layer is seen clearly.

151

 

 

that might explain why almost no literature discusses about the ring structure in

ferrofluids.

5.8 Conclusions

In this work we tried to study the early stages of the colloidal aggregation in
ferrofluid and the stability of the ring structure of these particles. We found that
the ring is the most stable structure for ferrofluid particles in low field and low
temperature as long as it has more than 4 particles. Based on simple arguments,
we derived an analytic expression for the energy of the ring and chain structure.
Two critical fields, 301 and Egg, are defined and their values from the conjugate-
gradient method was obtained. Performing the molecular dynamics simulation
including the quaternion parameters to represent the rotation of the dipoles, we
found the energy barriers, E‘, and the estimation of melting temperature TX, to
break a ring structure into lower symmetry ones as a function of applied field.
Our result indicates that in zero-field the ring structure should be observed in

carefully controlled experiments.

152

 

Bibliography

[Abro72]

[Alle84]

[Alle90]

[Ande62]

[Ande84]

[Arba90]

[Arca85]

Milton Abromowitz and Irene A. Stegan, Handbook of Mathematical

Functions, Dover Publications Inc., New York (1972).

Michael P. Allen, A molecular dynamics simulation study of octopoles

in the field of a planar surface, Mol. Phys. 52, 717 (1984).

M.P. Allen and DJ. Tildesley, Computer Simulations of Liquids, Ox-
ford University Press, New York (1990).

P.W. Anderson, Theory of flux creep in hard superconductors, Phys.

Rev. Lett. 9, 309 (1962).

H.C. Anderson, M.P. Allen, A. Bellemans, J. Board, J .H.R. Clarke,
M. Ferrario, J .M. Hail, S. Nosé,J.V. Opheusden and J .P. Ryckaert,
New molecular dynamics methods for various ensembles, Rapport

d’activité scientifique du CECAM, pp. 82-115 (1984)

S. Arbabi and M. Sahimi, Test of universality for three-dimensional

models of mechanical breakdown of disordered solids, Phys. Rev. B

41, 772 (1990).

L. de Arcangelis, S. Redner and H.J Herrmann, A Random fuse model

for breaking process, J. de Physique Lett. 46, L585 (1985).
153

 

[Ba1190] Pietro Ballone and Paolo Milani, Phys. Rev. B 42, 3201 (1990).

[Baro73] J. Barojas, D. Levesque and B. Quentrec, Simulation of diatomic

homonuclear liquids, Phys. Rev. A 7, 1092 (1973).

[Bea188] P.D. Beal and D.J. Srolovitz, Elastic fracture in random materials,
Phys. Rev. B 37, 5500 (1988).

[Bean62] C.P. Bean, Magnetization of hard superconductors, Phys. Rev. Lett.
8,250 (1962).

[Bean64] C.P. Bean, Rev. Mod. Phys. 36, 31 (1964).

[Beer82] N. Beer and D.G. Pettifor, Proc. NATO adv. study Inst. (Gent),
Plenum, New York (1982).

[Beng87] L. Benguigui, P. Ron and D.J. Bergman, Strain and Stress at the

fracture of percolative media, J. de Physique 48, 1547 (1987).

[Bil193] W.E. Billups and M.A. Ciufolini editors, Buckminsterfullerenes, VCH
Publishers, New York (1993).

[Bloc29] F. Bloch, Z. Phys. 52, 555 (1929).
[Bloc87] H. Bloch and J .P. Kelly, U.S. Patent No. 4687589 (1987).

[Bonn92] R.T. Bonnecaze and J .F. Brady, Dynamic simulation of an elec-

trorheological fluid, J. Chem. Phys. 96, 2183 (1992).

 

[Born91] I.C. van den Born, H.D. Hoekstra and J .'Th. de Hosson, Phys. Rev.

B 43, 3794 (1991).

154

[It

 

[Bowm87]

[Busm92]

[0.1.1094]

[Camp92]

[Camp93]

[Cast88]

[Chan91]

[Chau88]

[Chen92]

[ClarOO]

J.J. Bowman, T.B.A. Senior and P.L.E. Uslenghi, Electromagnetic
and Acoustic Scattering by Simple Shapes, Hemisphere Publishing
Corp., New York (1987).

H.C. Busmann, T. Lil], B. Reif and I.V. Hertel, Surf. Sci. 272, 146
(1992).

F. Leal Calderon, T. Stora, O. Mondain Monval, P. Poulin and J. Bi-
bette, Direct Measurement of Colloidal Forces, Phys. Rev. Lett. 72,

2959 (1994).

E.E.B. Campbell, A. Hielscher, R. Ehlich, V. Schyja and I.V. Her-
tel, Nuclear Physics Concepts in the Study Atomic Cluster Physics,
edited by R. Schmidt, H.C. Lutz and R. Dreizler, Lecture Notes in
Physics, Vol. 404, p. 185, Springer, Berlin (1992).

E.E.B. Campbell, V. Schyja, R. Ehlich and I.V. Hertel, Phys. Rev.
Lett. 70, 263 (1993).

E. Castillo, Extreme Value theory in Engineering, Academic Press,

San Diego, CA (1988).

S.K. Chan, S.G. Kim and RM. Duxbury, Elasticity and fracture of

tubular honeycombes containing random vacancies, Proc. CASME II

(1991).
P. Chaudhari et al., Phys. Rev. Lett. 60, 1653 (1988).
T.J. Chen, R.N. Zitter and R. Tao, Phys. Rev. Lett. 68, 2555 (1992).

David Clarke, private communication.

155

 

[II

 

 

[Che191]
[Clem87]

[0013989]

[Croc94]

[Cry089]

[Curl88]

[Curl91]

[Curt90]

[Curt91]
[Curt92]

[Dick93]

James R. Chelikowsky, Phys. Rev. Lett. 67, 2970 (1991).
J.R. Clem et al., Phys. Rev. B 35, 6637 (1987).

R.W. Coppard, L.A. Dissado, S.M. Rowland and R. Rakowski, Di-
electric breakdown of metal-loaded polyethylene, J. Phys.: Condens.
Matter 1, 3041 (1989).

John C. Crocker and David G. Grier, Microscopic Measurement of the
Pair Interaction Potential of Charge-Stabilized Colloid, Phys. Rev.

Lett. 73, 352 (1994).

Proceedings of the International Conference on “Critical Currents
in High Tc Superconductors”, Cryogenics, (29, 225-408 March 1989

supplement).
Robert F. Curl and Richard E. Smalley, Science 242, 1017 (1988).

Robert F. Curl and Richard E. Smalley, Scientific American, October
1991, p. 54.

W.A. Curtin and H. Scher, Brittle fracture of disordered materials,
J. Mater. Res. 5, 535 (1990).

W.A. Curtin and H. Scher, Phys. Rev. Lett. 67, 2457 (1980).
W.A. Curtin and H. Scher, Phys. Rev. B 45, 2620 (1992).

Akiva J. Dickstein, Shyamsunder Erramilli, Raymond E. Goldstein,
David P. Jackson and Stephen A. Langer, Labyrinthine Pattern For-

mation in Magnetic Fluids, Science. 261, 1012 (1993).

156

 

[Ii

 

[Dim088]

[Donn90]

[Dres88]

[Duxb86]

[Duxb87a]

[Duxb87b]

[Duxb90a]

[Duxb90b]

[Duxb91]

[Duxb94a]

 

D. Dimos et al., Phys. Rev. Lett. 61, 219 (1988).

J .B. Donnet and RC. Bansal, Carbon Fibers, Marcel Dekker Inc.,
2nd ed. (1990)

MS. Dresslhaus, G. Dresslhaus and K. Sugihara, I.L. Spain and
H.A. Goldberg, Graphite fibers and filaments, Springer Verlag, NY
(1988)

P.M. Duxbury, P.D. Beale and P.L. Leath, Phys. Rev. Lett. 57, 1052
(1986).

P.M. Duxbury, P.L. Leath and RD. Beale, Breakdown properties of
quenched random system — the random fuse network, Phys. Rev. B
36, 369 (1987).

P.M. Duxbury and P.L. Leath, J. Phys. A. 20, L411 (1987).

P.M. Duxbury and S.G. Kim, Scaling theory and simulations of frac-
ture in disordered media, Proc. CASME II (1990).

P.M. Duxbury, P.D. Beale, H. Bak and RA. Schroeder, Capacitance
and dielectric breakdown of metal loaded dielectrics, J. Phys. D. 23,
1546 (1990).

P.M. Duxbury and S.G. Kim, Scaling theory of elasticity and fracture
in disordered networks, Mat. Res. Symp. Proc. 207, 179 (1991).

P.M. Duxbury and P.L. Leath, Failure Probability and Average
Strength of Disordered Systems, Phys. Rev. Lett. 72, 2805 (1994).

157

 

[Duxb94b]

[Duze81]

[Ebbe92]

[Evan77a]

[Evan77b]

[Feng85]

[Feng87]

[Fili88]
[Foul89]

[Gala78]

[Ga1189a]

 

P.M. Duxbury, S.G. Kim and P.L. Leath, Size efiect and statistics of
fracture in random materials, Mat. Sci. and Eng. A176, 25 (1994).

T. Van Duzer and C.W. Turner, Principles of Superconductive De-

vices and Circuits, Elsevier, New York (1981)

T.W. Ebbesen, J. Tabuchi, K. Tanigaki, The mechanistics of fullerene
formation, Chem. Phys. Lett. 191, 336 (1992).

D.J. Evans, 0n the representation of orientation space, Mol. Phys.

34, 317 (1977).

D.J. Evans and S. Murad, Singularity-free algorithm for molecular

dynamics simulation of rigid polyatomics, Mol. Phys. 34, 327 (1977).

S. Feng, M.P. Thorpe and E. Garboczi, Efl'ective-medium theory of

percolation on central-force elastic networks, Phys. Rev. B 31, 276
(1985).

S. Feng, B.I. Halperin and RN. Sen, Transport properties of con-
tinuum systems near the percolation threshold, Phys. Rev. 35, 197
(1987).

RE. Filisko and W.E. Amstrong, U.S. Patent No. 4744914 (1988).
W.M.C. Foullers and R. Haydock, Phys. Rev. B 39, 12520 (1989).

J. Galambos, The Asymptotic Theory of Extreme Order Statistics,
Wiley, New York (1978).

G. Galli, R.M. Martin, R. Car and M. Parrinello, Phys. Rev. Lett.

62,555 (1989).
153

 

III

 

 

[c.0890]

[Garb88]

[Gear71]

[Gibs88]

[Gind92]

[Gold80]

[Good89]

[Grad80]

[Gumb58]

[Guyo87]

[11.1.90]

 

G. Galli, R.M. Martin, R. Car and M. Parrinello, Phys. Rev. Lett.
63, 988 (1989).

E.J. Garboczi, Efl'ective force constant for central force random net-

work, Phys. Rev. B 37, 318 (1988).

C.W. Gear, Numerical Initial Value Problems in Ordinary Differen-
tial Equations, Prentice-Hall, Englewood Cliffs, NJ (1971).

L.J. Gibson and M.F. Ashby, Cellular Solids, Pergamon Press, NY
(1988).

J .M. Ginder and L.D. Elie, Page 23, Electrorheological Fluids, edited
by R. Tao, World Scientific, Singapore (1992)

Herbert Goldstein, Classical Mechanics, 2nd Ed., Addison-Wesley,
Reading MA (1980).

L. Goodwin, A.J. Skinner and D.G. Pettifor, Europhys. Lett. 9, 701
(1989).

LS. Gradshteyn and I.M. Ryzhik, Table of integrals, series, and prod-
ucts, Academic Press, New York (1980).

E.J. Gumbel, Statistics of Extremes, Columbia Univ. Press, NY
(1958).

E. Guyon, S. Roux and D.J. Bergman, Critical Behavior of Elastic
failure thresholds in Percolation, J. de Physique 48, 903 (1987)

T.C. Halsey and W. Toor, Phys. Rev. Lett. 65, 2820 (1990).

159

 

[It

 

[Hans89]

[Har191]

[Hart77]

[Hash83]

[Hass89]

[Hayd80]

[Hein80]

[Held93]

[Herr89]

[Herr90]

[Hinr89]

 

A. Hansen, S. Roux and H.J. Herrmann, Rupture of central- orce

lattices, J. Phys. Paris 50, 733 (1989).
D.G. Harlow and S.L. Phoenix, J. Mech. Phys. Solids 39, 173 (1991).

EL. Hart, A Survey of the Literature on the size Efl'ect on Material
Strength, Tech. Rep. AFFDL-TR—77-11 (Air force Light Dynamics
Laboratory, Wright Patterson AFB, OH 45433), 1977.

Z. Hashin, Analysis of Composite Materials - A Survey, J. Appl.
Mech. 50, 481 (1983).

G.N Hassold and D.J. Srolovitz, Brittle fractureiin materials with

random defects Phys. Rev. B 39, 9273 (1989).

Roger Haydock, The Recursive Solution of the Schré'dinger Equation,
Solid State Phys. 35, 215 (1980).

Volker Heine, Electronic Structure from the Point of View of the Local

Atomic Environment, Solid State Phys. 35, 1 (1980).
G. von Helden, N.G. Gotts and M.T. Bowers, Nature 363, 60 (1993).

H.J. Herrmann, A. Hansen and S. Roux, Fracture of disordered elastic

lattices in two dimensions, Phys. Rev. B 39, 637 (1989).

H.J. Herrmann and S. Roux editors, Statistical Models for the Frac-
ture of disordered Materials, Elsevier Science Publishers, NY (1990)

EL. Hinrichsen, A. Hansen and S. Roux, A fracture growth model,
Europhys. Lett. 8, 1 (1989).

160

 

[It

 

[Hohe64]

[Hoov85]

[Hunt94]

01091]

[J ack75]

[J aya79]

[J oan74]

[J ohn90]

[J oos58]

[J und94]

utahnss]

[Keat66]

P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).

W.G. Hoover, Canonical dynamics: Equilibrium phase-space distri-

butions, Phys. Rev. A 31, 1695 (1985).

Joanna M. Hunter, James L. Fye, Eric J. Roskamp, and Martin
F. Jarrold, J. Phys. Chem. 98, 1810 (1994).

Sumio Iijima, Nature 354,, 56 (1991).

J .D. Jackson, Classical Electrodynamics, Second ed. Wiley, New York

(1975).

A. de S. Jayatilaka, Fracture of Engineering Brittle Materials, Ap-

plied Science, London (1979).
J. Joannopoulos and F. Yndurain, Phys. Rev. B 10, 5164 (1974).

Robert Johnson, Gerard Meijer and Donald S. Bethune, 060 Has
Icosahedral Symmetry, J. Am. Ceram. Soc. 112, 8983 (1990).

G. J oos, Theoretical Physics, 3rd ed., p. 218, Haffner Publishing Com-
pany, New York, N.Y. (1958).

P. Jund, S.G. Kim, D. Tomének and J. Hetherington, Stability of

Magnetic Dipole Structures in Ferrofluids, preprint (1994).

B. Kahng, G.G. Bartrouni, S. Redner, L. de Arcangelis and H.J. Her-
rmann, Electric breakdown in a fuse network with continuously dis-

tributed breaking strengths, Phys. Rev. B 37, 7625 (1988).

P.N. Keating, Phys. Rev. 145, 637 (1966).
161

 

[Kell80]

[Ke1181]

[Kell86]

[Kern92]

[KimE93]

[KimS91]

[KimS 94a]

[KimS 94b]

[Kim894c]

[KimY62]

[Kirk39]

M.J Kelly, Application of the Recursion Method to the Electronic
Structure from an Atomic Point of View, Solid State Phys. 35, 295
(1980).

B.T. Kelly, Physics of Graphite, Applied Science Publishers (1982).

A. Kelly and NH. McMillan, Strong Solids, Clarendon Press, Oxford
(1986).

R. Kerner, K.A. Penson and K.H. Bennemann, Europhys. Lett. 19,

363 (1992).

Eunja Kim, Young Hee Lee and Jae Young Lee, Fragmentation of

060 and C70 clusters, Phys. Rev. B 48, 18230 (1993).

S.G. Kim and P.M. Duxbury, Cracks and critical current, J. Appl.

Phys. 70, 3164 (1991).

Seong Gon Kim and David Tomanek, Melting the fullerenes: A

molecular dynamics study, Phys. Rev. Lett. 72, 2418 (1994).

S.G. Kim, W. Zhong and D. Tomének, Synthesizing Carbon

Fullerenes, preprint (1994).

S.G. Kim, P. Jund, D. Tomanek and J. Hetherington, Structure For-

mation in Ferrofluids, in preparation (1994).

Y.B. Kim, C.F. Hempstead and A. Strnad, Phys. Rev. Lett. 9, 306

(1962).
J.G. Kirkwood, J. Chem. Phys. 7, 506 (1939).

162

 

[Kitt86]

[Klin89]

[Kohn65]

[Kraj82]

[Krat90]

[Krot85]

[Kusa90]

[Lang92]

[Leat94]

[Len092]

[ng3)

C. Kittel, Introduction to Solid State Physics, sixth edition, John
Wiley, New York (1986).

D.J. Klingenberg, Frank van Swol and C.F. Zukoski, Dynamic sim-
ulation of electrorheological suspensions, J. Chem. Phys. 91, 7888
(1989).

W. Kohn and L.J. Sham, Phys. Rev. 140, A1133 (1965).

D. Krajcinovic and M.A.G. Silva, Statistical Aspects of the Continu-

ous Damage Theory, Int. J. Solids Structures 18, 551 (1982).

W. Kritschmer, L.D. Lamb, K. Fostiropoulos, and DR Huffman,

Nature 347, 354 (1990).

H.W. Kroto, J.R. Heath, S.C. O’Brien, R.F. Curl, and RE. Smalley,
Nature 318, 162 (1958).

RC. Kusalik, Computer simulation results for the dielectric proper-

ties of a highly polar fluid, J. Chem. Phys. 95, 3520 (1990).

Stephen A. Langer, Raymond E. Goldstein and David P. Jackson, Dy-
namics of labyrinthine pattern formation in Magnetic Fluids, Phys.

Rev. A 46, 4894 (1992).

P.L. Leath and P.M. Duxbury, preprint.

Thomas Lenosky, Xavier Gonze, Michael Teter and Veit Elser, Nature
355, 333 (1992).

T. Lill, F. Lacher, H.-G. Busmann and I.V. Hertel, Phys. Rev. Lett.

70, 3383 (1993).
163

 

[LiY887]

[LiY888]

[LiYS89]

[LiXP93]

[Loui87]

[Luch87]

[Lund83]

[Mann88]
[Mann89]

[Mand84]

[Mart92]

 

Y.S. Li and P.M. Duxbury, Size and location of the largest current in

a random resistor network, Phys. Rev. B 36, 5411 (1887).

Y.S. Li and P.M. Duxbury, Crack arrest by residual bonding in spring
and resistor networks, Phys. Rev. B 38, 9257 (1988).

Y.S. Li and P.M. Duxbury, Phys. Rev. B 40, 4889 (1989).

X.-P. Li, R.W. Nunes and D. Vanderbilt, Phys. Rev. B 47, 10891
(1993).

E. Louis and F. Guinea, The fractal nature of fracture, Europhys.
Lett. 3, 871 (1987).

M.U. Luchini and C.M.M. Nex, A new procedure for appending ter-

minators in the recursion method, J. Phys. C: Solid State Phys. 20,

3125 (1987).

S. Lundqvist and N.H. March, Theory of the Inhomogeneous Electron

Gas, Plenum (1983).
J. Mannhart et al., Phys. Rev. Lett. 61, 2476 (1988).
J. Mannhart et al. Science. 245, 839 (1989).

BB. Mandelbrot, D.E. Passoja and A.J. Paullay, Fractal character

of fracture surfaces of metals, Nature 308, 721 (1984).

J .E. Martin, J. Odinek and T.C. Halsey, Phys. Rev. Lett. 69, 1524

(1992).

164

 

(‘3

 

[McLa46]

[McLa47]

[Meak88]

[Meak91]

[Mech89]

[Mors53]

[Mull49a]
[Mull49b]

[Nels78] .

[Ne1388]

[Nick88]

N.W. McLachlan, Mathieu functions and their classification, J. Math.
Phys. 25, 209 (1946).

N .W. McLachlan, Theory and Application of Mathieu Functions,
Clarendon Press, Oxford, England (1947).

P. Meakin, Simple models for crack growth, Crystal prop. and prep.
17 and 18, 1, (1988).

P. Meakin, Science 252, 226, (1991).

J .J . Mecholsky, D.E. Passoja and K.S. Feinberg-Ringel, Quantitative
analysis of brittle fracture surfaces using fractals geometry, J. Am.

Ceram. Soc. 72, 60 (1989).

Philip M. Morse and Herman Feshbach, Methods of Theoretical
Physics, Part I and Part II, McGraw-Hill Book Company, Inc., New
York (1953).

R.S. Mulliken, J. Chem. Phys. 46, 467 (1949).
R.S. Mulliken, J. Chem. Phys. 46, 675 (1949).

David R. Nelson, Study of melting in two dimensions, Phys. Rev. B

18, 2318 (1878).
DR. Nelson, Phys. Rev. Lett. 60, 1973 (1988).

G.A. Nicklasson, A. Torebring, C. Larsson, G.G. Granqvist and
T. Farestam, Fractal Dimension of Gas-Evaporated C’o Aggregates:

Role of Magnetic Coupling, Phys. Rev. Lett. 60, 1735 (1988).

165

 

’11

 

[Nose84]

[Orla91]

[Page52]

[Phoe83]
[Phoe92]

[Piggsm

[Pres86]

[RayP85]

[Reif65a]

[Reif65b]

[Rose85]

S. Nosé, Mol. Phys. 52, 255 (1984).

Terry P. Orlando and Kevin A. Delin, Foundation of Applied Super-
conductivity, Addison Wesley, Reading MA (1991).

L. Page, Introduction to Theoretical Physics, 3rd ed., p. 286, D. Van

Nostrand Company, Princeton, NJ. (1952).
S.L. Phoenix and L.J. Tierney, Eng. Frac. Mech. 18, 193 (1983).
S.L. Phoenix and R.L. Smith, preprint (1988).

M.R. Piggot, Load Bearing Composites, Pergamon Press, Oxford
(1980).

William H. Press, Brian P. F lannery, Saul A. Teukolsky and William
T. Vetterling, Numerical Recipes: The Art of Scientific Computing,
Cambridge University Press, Cambridge (1986).

P. Ray and B.K.Chakrabarti, A microscopic approach to the statisti-
cal fracture analysis of disordered brittle solids, Solid State Commu-

nications 53, 477 (1985).

Federick Reif, Fundamentals of Statistical and Thermal Physics, Page
248, McGraw-Hill, New York, NY (1965).

Federick Reif, Fundamentals of Statistical and Thermal Physics, Page
560, McGraw-Hill, New York, NY (1965).

RE. Rosensweig, Ferrohydrodynamics, Cambridge University Press,
Cambridge (1985).

166

 

[Root51]

[Root60]

[Roux88]

[Sahi86]

[Sahi91]
[3:11.193]

[Seif92]

[Sier00]

[Sier86]

[SihG73]

[Slat54]

[Smit80]

[Sr0188]

C.C.J. Roothaan, Rev. Mod. Phys. 23, 69 (1951).
C.C.J. Roothaan, Rev. Mod. Phys. 32, 179 (1960).

S. Roux, A. Hansen, H.J. Herrmann and E. Guyon, Rupture of Het-
erogeneous Media in the limit of infinite disorder, J. Stat. Phys. 52,
237 (1988).

M. Sahimi and J .D. Goddard, Elastic percolation models for cohesive

failure in heterogeneous systems, Phys. Rev. B 33, 7848 (1986).
M. Sahimi and S. Arbabi, Mat. Res. Symp. Proc. 207, 201 (1991).
M. Sahimi and S. Arbabi, Phys. Rev. B 47, 713 (1993).

G. Seifert and R. Schmidt, p. 261, Clusters and Fullerenes, edited by
V. Kumar, T.P. Martin and E. Tosatti, World Scientific, Singapore

(1992).
K. Sieradzki, private communication.

K. Sieradzki and R. Li, Fracture behavior of a solid with random

porosity, Phys. Rev. Lett. 56, 2509 (1986).

G.C. Sih, Handbook of Stress Intensity Factors, Institute of Fracture

and Solid Mechanics, Lehigh University, Bethlehem, PA (1973).

J .C. Slater and G.F. Koster, Simplified LCAO Method for the Peri-

odic Potential Problem, Phys. Rev. 94, 1498 (1954).
R.L. Smith, Proc. R. Soc. London A 372, 539 (1980).

D.J. Srolovitz and RD. Beale, J. Am. Ceram. Soc., 71, 362 (1988).
167

 

[Stev94]

[Tada73]

[Tak386]

[TaoR92]

[TaoR94]

[Teje94]

[Term86]

[Term87a]
[Term87b]

[Term88]

Mark J. Stevens and Gary S. Grest, Coexistence in Dipolar Fluids in
a Field, Phys. Rev. Lett. 72, 3686 (1994).

H. Tada, P.C. Paris and G.R. Irwin, The Stress Analysis of Cracks
Handbook, Del Research Corp., Hellertown, PA (1973).

H. Takayasu, Pattern formation of dentritic fractals in fracture and
electric breakdown, Proceedings of the symposium on fractals in

physics, 181 (1986).

R. Tao, editor, Electrorheological Fluids, World Scientific, Singapore

(1992)

R. Tao and Qi Jiang, Simulation of Structure Formation in an Elec-

trorheological Fluid, Phys. Rev. Lett. 73, 205 (1994).

C.F. Tejero, a. Daanoun, H.N.W.Lekkerkerker and M. Baus, Phase
Diagrams of “Simple” Fluids with Extreme Pair Potentials, Phys.

Rev. Lett. 73,752 (1994).

Y. Termonia and P. Meakin, Formation of fractal cracks in a kinetic

fracture model, Nature 320, 429 (1986).
Y. Termonia and P. Smith, Macromolecules, 20, 835 (1987).
Y. Termonia, J. Mater. Sci. 22, 1733 (1987).

Y. Termonia and P. Smith, in A.E. Zachariades and R.S. Porter (ed-
itors), “High Modulus Polymers Approaches to Design and Develop-

ment”, Dekker, New York (1988).

168

 

’11

 

[Ters88]
[Thor83]

[Thor8 7]

[Ti1186]

[Tink65]

[Toma86]

[Toma91]

[Turc85]

[Unde86]

[Vand92]

[wmgz]

J. Tersofi', Phys. Rev. Lett. 61, 2879 (1988).
M.F. Thorpe, J. Non-Cryst. Solids 57, 355 (1983)

M.F. Thorpe, Rigidity Percolation in Glassy Structures, J. Non-
Cryst. Solids 76, 109 (1987)

DR. Tilley and J. Tilley, Superfluidity and Superconductivity, Adam
Hilger Ltd., Boston (1986).

M. Tinkham, Superconductivity, Gordon and Breach, New York
(1965)

David Tomének and Michael A. Schluter, Calculation of Magic Num-
bers and the Stability of Small Si Clusters, Phys. Rev. Lett. 56, 1055

(1986).

David Tomanek and Michael A. Schluter, Growth Regimes of Carbon
Clusters, Phys. Rev. Lett. 67, 2331 (1991).

D.L Turcotte, R.F. Smalley Jr. and Sara A. Solla, Collapse of loaded

fractal trees, Nature 313, 671 (1985).

E.E. Underwood and K. Banerji, Fractals in fractography, Mat. Sci.

and Eng. 80, 1 (1986).
David Vanderbilt and J. Tersoff, Phys. Rev. Lett. 61, 2879 (1988).

T. Wakabayashi and Y. Achiba, A model. for the 060 and 070 growth

mechanism, Chem. Phys. Lett. 190, 465 (1992).

169

 

[Waka93]

[Wang94]

[WeasQO]

[Weib51]
[WeiD92a]

[WeiD92b]

[Wei593]

[Welt89]

[XiaW89]

[XuCH92]

T. Wakabayashi,H.Shiromaru, K. Kikuchi and Y. Achiba, A selective

isomer growth offullerenes, Chem. Phys. Lett. 201, 470 (1993).

Hao Wang, Yun Zhu, C. Boyd, Weili Luo, A. Cebers and
RE. Rosensweig, Periodic Branched structures in a Phase-Seperated

Magnetic Colloid, Phys. Rev. Lett. 72, 1929 (1994).

Robert C. Weast, editor in-chief, CRC Handbook of Chemistry and
Physics, Page B-10, 62th edition, CRC Press, Boca Raton, Florida
(1990).

W. Weibull, J. Appl. Mech. 18, 293 (1951).
Dongquing Wei and G.N. Patey, Phys. Rev. Lett. 68, 2043 (1992).

Dongquing Wei and G.N. Patey, Ferroelectric liquid-crystal and solid
phases formed by strongly interacting dipolar soft spheres, Phys. Rev.
A 46, 7783 (1992).

J .J . Weis and D. Levesque, Chain Formation in Low Density Dipo-
lar Hard Spheres: A Monte Carlo Study, Phys. Rev. Lett. 71, 2729

(1993).

William Weltner, Jr. and Richard J. Van, Chem. Rev. 89, 1713
(1989).

W. Xia and P.L. Leath, Phys. Rev. Lett. 63, 1428 (1989).

C. H. Xu, C. Z. Wang, C. T. Chan and K. M. Ho, A transferable
tight-binding potential for carbon, J. Phys.: Condens. Matter 4, 6047

(1992).
170

 

' [XuCH94] Chunhui Xu and Gustavo E. Scuseria, Tight-Binding Molecular Dy-
namics Simulations of Fullerene Annealing and Fragmentation, Phys.

Rev. Lett. 72, 669 (1994).

[Yann91] C.S Yannoni, P.P Bernier, D.S. Bethune, G. Meijer and J .R. Salem,
J. Am. Ceram. Soc. 113, 3190 (1991).

[Zab086] J.G.Zabolotzky, D.J. Bergman and D. Stauffer, J. Stat. Phys. 44,
211 (1986).

[Zhon93] W. Zhong, D. Tomének and George F. Bertsch, Total energy calcula-
tions for extremely large clusters: The recursive approach, Solid State

Communications 86, 607 (1993).

171

 

l
88

will“

LLQHOZG

l l
1

V
T.
N
U
E
T
on
T
S
N
AH
G
T.
H
C
T.
H