3.,3. J

o.
v
v.

 

MS U is an Affirmative Action/Equal Opportunity Institution

MICHI AN STATE UNIVERSI

lull/ll“ x will

my ll

ll 1/ HM“

This is to certify that the

dissertation entitled

STRUCTURE AND DYNAMICS OF
LAYERED SYSTEMS

presented by

Hyangsuk Seong

has been accepted towards fulfillment
of the requirements for

Pk. D dew“, 7>H YSICS

 

JDW

Major professor

Dateidpmt 21, 1995

 

0-12771

 

AH —.——*——_ r - ‘7‘

i, —\
LIBRARY
Michigan State
University

 

 

L

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

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

_|L'

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU Is An Affirmative ActionEqual Opportunity Institution
encirchma-pd

STRUCTURE AND DYNAMICS OF
LAYERED SYSTEMS

By

Hyangsuk Seong

A DISSERTATION

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

DOCTOR OF PHILOSOPHY

Department of Physics and Astronomy

1993

‘

ABSTRACT

Structure and Dynamics of Layered Systems

BY

Hyangsuk Seong

This dissertation addresses some of the current problems dealing with the struc-
ture and dynamics of layered solids. Three different types of physical systems have
been studied. (1) 2—dimensional (2D) repulsive screened Coulomb systems in the
presence of corrugation modelling a large class of intercalation compounds, (2) ran-
domly intercalated layered solids modelling ternary intercalation systems, and (3)

layered superlattices with interfacial disorder.

For (1), using molecular dynamics simulations we have investigated structural,
thermodynamic, and dynamic properties as a function of the intercalant density and
the strength of the corrugation. For physical parameters appropriate for RbCu, we
explain the observed low-T structure of the Rb intercalants, the continuous nature
of the melting transition, and the wave number and frequency dependence of the
dynamic structure factor observed in neutron scattering measurements. In addition,
we have studied in detail how the dynamics of a 2D fluid evolves from homogeneous
(corrugation free) to lattice-fluid limit by changing the strength of the corrugation.
We find that corrugation enhances the possibility of observing two—phonon peak
in the solid phase and these peaks persist even in the liquid phase, a remarkable

property of a lattice—fluid.

For (2), we have studied the combined effects of local anharmonicity and the

transverse layer rigidity on the gallery expansion. Simulation results are compared

with those obtained within effective medium approximation and a simple and yet
successful model called the Catchment Area Model. While in the harmonic limit
Vegard’s law is obtained when the two intercalants have the same compressibility,
we find that inclusion of anharmonicity gives rise to deviations from Vegard’s law

which increase with the host layer rigidity and the degree of anharmonicity.

For (3), the effect of interfacial disorder on the structure of coherent superlat-
tices has been studied. For a given profile of the interlayer diffusion, various physical
quantities such as average bond lengths, specific bond lengths, and change in the
coherent strain energy with the degree of interfacial disorder have been calculated.
We show that depending on local geometry, interfacial disorder can in fact be ener-

getically favorable.

To My Parents

ACKNOWLEDGMENTS

I deeply acknowledge my advisor, Professor S.D. Mahanti, for his continual and
patient guidance throughout all stages of this work. Without his understanding,

support, and influence on my graduate study this dissertation could not have been

finished.

I thank Professor M.F. Thorpe for his advice and collaboration on the problem
of disordered interfaces in coherent multilayer systems and the intercalated bilayer
problems. I also thank Professor J .S. Chung for collaboration on the problem of
gallery expansion of randomly intercalated anharmonic bilayers. I would like to
thank Drs. Surajit Sen and Tahir Cagin for their advice and collaboration on the
structure of 2-dimensional intercalation compounds. I specially thank Dr. Tahir

Cagin for giving me his MD code and for his computer programming skill.

I am grateful to Professors M.F. Thorpe, T.A. Kaplan, C. Foiles, and D. Stump
for their service on my Ph.D. guidance committee. I would like to thank Professor

J .J. Kovacs, Ms. J. King, and Ms. S. Holland for their assistance.

My special thanks go to my parents for their support and encouragement. I
also thank my dear friends: Junghwa Park Bang, Hyunju Chang, Keng Ching Lin
and others in this department for their friendship and refreshing conversation. I

acknowledge Seonggon Kim for his help in solving all kinds of problems encountered

in VAX.

Financial support from Michigan State University, the center for Fundamantal

Material research and the National Science Foundation are gratefully acknowledged.

iv

Contents

1 Introduction 1
2 Structure of Disorderd Interfaces in Coherent Multilayers 9
2.1 Introduction ................................ 9
2.2 The Model ................................. 10
2.3 Analytic Solutions and Numerical Simulations ............. 13
2.3.1 Average Bond Lengths ...................... 15

2.3.2 Specific Bond Lengths ...................... 20

2.3.3 Fluctuations and Strain Energy ................. 23

2.4 Summary and Conclusion ........................ 26

3 Gallery Expansion of Randomly Intercalated Anharmonic Bilayers 31

3.1

3.2

3.3

3.4

3.5

Introduction ................................ 31
The Model and Limiting Cases ..................... 35
Single Impurity with Anharmonic Cubic Potential ........... 40
Effective Medium Theory (EMT) .................... 42
Numerical Simulations and Comparison with EMT .......... 48

3.6 Comparison of Anharmonic Potential Model with Catchment Area

3.7 Conclusion ................................. 60

Structure of 2-Dimensional (2D) Repulsive Screened Coulomb Sys-

tems/ Model Intercalation Compounds 63
4.1 Introduction ................................ 63
4.2 Physical Systems ............................. 64
4.3 Potential Model .............................. 68
4.4 Details of Molecular Dynamics Simulations ............... 69
4.5 Structural Properties ........................... 70
4.6 Periodic Domain Wall Model ....................... 76
4.7 Summary and Conclusion ........................ 81

Melting of a Repulsive Screened Coulomb System in 2D - Effect of

Corrugation 85
5.1 Introduction ................................ 85
5.2 Physical Systems ............................. 88
5.3 Thermodynamic PrOperties ....................... 89
5.3.1 Energy vs. Temperature ..................... 90
5.3.2 Bond Orientational Order Parameter .............. 92
5.3.3 Angular Susceptibility ...................... 94
5.3.4 Topological Defects ........................ 97

vi

5.3.5 Translational Diffusion Constants ................ 99

5.4 Summary and Discussion ........................ 103

Dynamics of 2D Repulsive Screened Coulomb Fluids on a Corru-

gated Surface 109
6.1 Introduction ................................ 109
6.2 Time Correlation Functions in Real Space ............... 111
6.2.1 G.(r,t) for Typical Systems ................... 112
6.2.2 Self-Diffusion Constant ...................... 114
6.3 Intermediate Scattering Function .................... 121
6.4 Incoherent Dynamic Structure Factor .................. 128
6.4.1 Classical Harmonic Oscillator .................. 131
6.4.2 Corrugation-free System ..................... 134
6.4.3 Effect of a Strong Corrugation Potential ............ 138
6.4.4 Evolution of S,(q, u) with Strength of Corrugation ....... 141
6.5 Coherent Correlation Functions ..................... 144
6.6 Summary and Discussion ......................... 157
Multiphonon Response in 2D Corrugated Systems 161
7.1 Introduction ................................ 161
7.2 Lattice Dynamics ............................. 166

7.3 Phonon Density of States through Velocity Auto—correlation Function 170
7.4 Phonon Density of States through Incoherent Dynamic Structure F actor174

vii

7.5 Summary and Conclusion ........................

A Average length for A at O
A.1 6 - bonds (same layer) . .' ........................

A.2 11(1)) - bonds (interlayer) .........................

B EMT in Harmonic Limit

C Mean Square Displacement in d-dimension

viii

182

185

185

186

188

189

List of Tables

4.1 Epitaxy angles and primary peak positions of S (k) .......... 76

6.1 Barrier height and diffusion constant for various systems ....... 117

ix

List of Figures

1.1

2.1

2.2

2.3

2.4

2.5

2.6

2.7

3.1

3.2

3.3

3.4

3.5

3.6

A schematic phase diagram in the corrugation strength vs. temperature 3

Heterostructures with interfacial disorder ................ 11
Schematic picture of coherent superlattices ............... 16
Average layer spacings .......................... 18
Average bond length vs. 2:5 near the interface layers .......... 19
Specific bond lengths of the interface in the 6—direction ........ 22
Specific bond lengths near the interfaces in the u-direction ...... 24
Strain energy as a function of interfacial disorder ........... 27
Schematic pictures of bilayers with different layer rigidities ...... 34
Typical average spacings for different interaction potentials ...... 39

Comparison of the numerical simulations with a single impurity results 43
Average gallery spacing of a bilayer system ............... 50

Comparison of the numerical simulations with the results of EMT for

different layer rigidities .......................... 51

Comparison of the numerical simulations with the EMT results for

different strengths of the anharmonicity ................. 53

3.7

3.8

3.9

4.1

4.2

4.3

4.4

4.5

4.6

5.1

5.2

5.3

5.4

5.5

5.6

5.7

6.1

Comparison of the simulation results obtained with the LJ potential
model with the Catchment Area Model ................. 55

Comparison of the simulation results obtained with the LJ potential

model with the modified Catchment Area Model I ........... 58

Comparison of the simulation results obtained with the LJ potential

model with the modified Catchment Area Model 11 .......... 59
Structure of graphite and staging sequence ............... 65
Pair correlation function g(r) for R120“ ................. 72
Solid structure of RbC'24 and Rng4,57 .................. 73
Static structure factor 3 (k) for R5024 at 3K .............. 75
Static structure factor for Rb024 at 250K ................ 77
Two epitaxy angles vs. various domain sizes .............. 80
Potential energy vs temperature ..................... 91

Bond orientational order parameter vs temperature for systems I and II 93
Angular susceptibility vs temperature for systems I and II ...... 95

Temperature dependence of (a) bond orientational order parameter

for systems II and III and (b) corresponding angular susceptibility . . 98

Topological defects in solid and liquid states .............. 100
Number of disclinations vs temperature ................. 101
Diffusion constant vs temperature .................... 102

Trajectories of tagged particles for various strengths of corrugation . . 115

xi

6.2

6.3

6.4

6.5

6.6

6.7

6.8

6.9

6.10

6.11

6.12

6.13

6.14

6.15

6.16

6.17

6.18

Time dependent diffusion coefficient near 250K ............ 116

VAF for various values of 2K near 250K ................ 119
VAF for various temperatures for 2K = 0.9 ............... 122
The self part of the intermediate scattering functions at q = 0.4;?1 . 124

The self part of the intermediate scattering function at q = 1.2:?1 . . 125
The self part of the intermediate scattering function at q = 4A“1 . . . 126

Spectral weight of the incoherent dynamic structure factor of a har-

monic oscillator at u = 0 for various values of 2K ........... 132

Spectral weight of the incoherent dynamic structure factor of a har-

monic oscillator at V = uo for various values of 2K ........... 133
Incoherent dynamic structure factor vs. frequency for a system with-
out corrugation potential at 254K .................... 135
Diffusion coefficients from S.(q, u) as functions of q .......... 137
Incoherent dynamic structure factor in the liquid state vs. frequency

for a strong corrugated potential (2K = 0.9) .............. 139

Incoherent dynamic structure factor vs. frequency for various strengths

of the corrugation near 250K for q = 424‘1 ............... 142
Intermediate scattering function at q = 0.4?1‘1 near 250K ....... 146
Intermediate scattering function at q = 1.22?1 near 250K ....... 147
Intermediate scattering function at q = 4A" near 250K ........ 148
Dynamic structure factor at q = 0AA" ........ i ......... 150
Dynamic Structure factor at q = 1.2:?1 ................ 151

xii

6.19 Dynamic Structure factor at q = 4104'1 ................. 152

6.20 Longitudinal current correlation functions for 2K = 0 ......... 154
6.21 Longitudinal current correlation functions for 2K = 0.9 ........ 155
6.22 Dispersion curve of Brillouin peaks ................... 156
7.1 Intercalant in-plane mode phonon density of states from Kamitaka-

hara and Zabel for R502. ......................... 163
7.2 Unit cells of systems I and III ...................... 168
7.3 Dispersion curves for systems I and III ................. 169
7.4 Phonon density of states through lattice dynamics for the system III . 171
7.5 Fourier transform of velocity auto-correlation function for the systems

II and III ................................. 172
7.6 A(q, u) at q = 4A“ for the system II .................. 176
7.7 A(q, u) at q = 0.4;!”1 for the system II ................. 178
7.8 A(q, u)’sifor a harmonic oscillator .................... 180
7.9 Multiphonon response of the systems I and II ............. 181

xiii

Chapter 1

Introduction

A large number of solids occurring in nature or synthesized in a laboratory exhibit
a high degree of anisotropy in their physical properties. In some of these solids,
the interatomic forces within a layer are usually much stronger than those between
the layers. Such solids are characterized as layered solids, typical examples are
graphite, layered chalcogenides, and silicate clays. In other types of layered solids,
although interatomic forces are not so anisotropic, certain physical properties are
highly so. Examples are high-Tc oxide superconductors and artificial superlattices.
Many of these layered solids provide a natural arena to study two dimensional (2D)

or quasi-2D physical phenomena.

This dissertation focuses on some of the physical properties of layered solids; (1)
the structure, melting, and dynamics of 2D systems in the presence of substrate
corrugation modelling the properties of binary intercalation systems (chapters 4-7),
(2) gallery spacing of ternary intercalation systems (chapter 3), and (3) the struc-
ture of layered superlattices in the presence of interfacial disorder (chapter 2). The
layered solids with weak interlayer interaction can be made to imbibe guest molec-
ular, ionic and/or atomic species which are sandwiched between the host layers.

This process is called intercalation and manifests itself most dramatically as a gross

one-dimensional expansion of the pristine host in a direction normal to the layer
planes (denoted commonly as the c-axis). For (1) and (2), the systems we have in
mind are binary and ternary intercalated layered solids, the physical properties of

the intercalants being of interest in (1).

In binary intercalation compounds, where there is only one type of intercalant
inside the gallery, the intercalants form a two dimensional sheet and their structure
is determined by a competition between the interaction between the intercalants and
the intercalant-host interaction. Depending on the planar density of the intercalants
and the strength of the substrate corrugation, one expects to see incommensurate
structures characterized by domains and domain walls of various size, or simple
commensurate structures. Thermodynamic properties, in particular melting of these
2D solids, and their dynamic properties both in solid and liquid states are areas of

great current interest. We address these questions in the present thesis.

In Figure 1.1, we show a generic phase diagram of a 2D system in the presence
of corrugation, K being a measure of its strength. Of course the nature of the in-
terparticle potential whether repulsive screened Coulomb (charged intercalants or
ionic overlayers) or Lennard Jones types (rare gas atoms on graphite and similar
systems), and the density vis-a-vis the period of corrugation play extremely impor-
tant roles in the detailed nature of the phase diagram. When K = 0, one has the
problem of melting of a homogeneous 2D solid. In the limit of extremely large K,
the particles are localized at the minima of the corrugation potential thereby giving
a “lattice-gas” system whose thermodynamic properties are quite different from the
K = 0 case. We will focus on two aspects of this phase diagram for a repulsive
screened Coulomb system; these are shown as two lines in the figure. First, for

a given strength of the corrugation, we want to see how the system melts from a

 

 

g s
n
=3 3
3° A ‘
s i
-. a
i a
h n
on

Solid Liquid

 

 

 

Temperature ('1')

Figure 1.1: A schematic phase diagram in the corrugation strength (K) vs. temper-
ature (T) plane is given for a constant density. The solid curve indicates transition
from a solid to a liquid. The study along a dotted line denoted as A gives the T
dependence of melting for a fixed strength of the corrugation. The study along the
dashed line, B, tells how the liquid dynamics changes with K at a given T.

solid to a liquid and how this differs from the K = 0 case (A). Second, for a given
temperature, we study how the structure and dynamics of a 2D fluid evolve from a

homogeneous fluid to a lattice-fluid (B).

In ternary intercalation compounds, one has a mixture of two types of inter-
calants inside the gallery. In general, structural, thermodynamic, and dynamic
properties of these systems are quite complex. One problem that has been of great
interest in recent years is the overall gallery structure of ternary intercalation com-
pounds of the type A1-,B,L, where A and B are two types of intercalants of different
size and L is the host. One studies the average c-axis spacing as a function of a:

for a random configuration of the intercalants and tries to understand the general

conditions under which one sees a linear a: dependence (Vegard’s law) and devia-
tions from this law. In; addition, there have been several interesting studies on the
ground state structure of 2D random alloys of two different types of atoms with
different sizes to see how the system looses its long range order and algebraic de-
cay in correlation as either a function of size mismatch or alloy concentration, the
so called size-mismatch melting. Here, we are concerned only with the problem of
gallery expansion with specific purpose of identifying the role of anharmonicity and

its effects on the Vegard’s law.

In coherent superlattices with interfacial disorder, another layered system of
great current interest, it is believed that disorder near the interfaces is very im-
portant and may profoundly affect the physical properties of these systems. In
this thesis, our primary interest is to study the effect of interfacial disorder on the

structures of these superlattice systems.

As a first step towards understanding the interfacial structure of coherent multi-
layer systems, we use harmonic potentials for the interaction of atoms in multilayers
(chapter 2). Next, to study the gallery structure of randomly intercalated ternary
systems (chapter 3), anharmonic potentials are introduced for the interaction be-
tween intercalants and host layers while keeping harmonic potentials for interaction
between the host atoms themselves. In these two chapters, we address the T=0
structure, the disorder configuration being frozen. From chapters 4 to 7, we use
anharmonic potentials for both intercalant (adsorbate)-host (substrate) interactions
and intercalant-intercalant interactions to study the structure and dynamics as a
function of temperature. In the following, a brief description of the physical sys-

tems and their properties which we have investigated is presented.

Chapter 2 is the study of the effect of interfacial disorder on the structure of

coherent superlattices consisting of 1 layers of A atoms and m layers of B atoms. For
a given profile of the interlayer diffusion, i.e the concentration of A atoms in B rich
layers and vice versa, we have been able to obtain several exact results for different
structural parameters characterizing the disordered systems in the harmonic limit if
we assume that the spring constants associated with the straining of AA, AB, and
BB bonds are the same. Various physical quantities such as (1) the average bond
lengths both parallel and perpendicular to the growth direction, (2) specific bond
lengths such as AA, AB, and BB both near the interface and away from it, and
(3) the change in the coherent strain energy with degree of interfacial disorder have
been calculated. Our theoretical results can be applied to the multilayers of Ge/Si,
GaAs/GaInAs, and GaAs/GaSbAs.

Chapter 3 concerns with the gallery expansion of randomly intercalated anhar-
monic systems. In these systems one usually has a random alloy of two types of
intercalants inside the gallery of host layers, denoted as A1-,B,L, and the gallery
expansion is monitored as a function of 2:. These systems are referred to as ternary
intercalation compounds. Harmonic models have been extensively studied in bilay-
ers, multilayers, and in 2D and 3D solid solutions. General conditions for the validity
'of Vegard’s law within harmonic theory have been put forward. Whereas some 3D
solid solutions show Vegard’s law behavior, randomly intercalated layered solids
show significant deviations from Vegard’s law. It is believed that anharmonicity ef-
fects in layered intercalation compounds are significant. Therefore, we have studied
the combined effects of local anharmonicity and the transverse layer rigidity on the
gallery eXpansion in ternary intercalated systems, but for the simplicity we have
chosen a bilayer model. Unlike the harmonic model, the problem cannot be solved

exactly when anharmonicity is present. Numerical simulation results are compared

with analytic calculations within an effective medium approximation and a simple

but quite successful “Catchment Area Model”.

As discussed above, chapters 2 and 3 are primarily studies of the structure of
partially disordered or disordered solids at T = 0. For a given random configuration
of the intercalants their positions are allowed to relax and the role of disorder on
the structure is the main concern. On the other hand, chapters 4 through 7 involve
a study of structural, melting, and dynamic properties of a multicomponent system
where one of the components is frozen due to its own structural rigidity and the other
forms a 2D physical system and shows many interesting phenomena. For example
in layered intercalated compounds, the host layer is assumed to be rigid and the
intercalants, in addition to a screened Coulomb repulsion, feel the static corrugation
potential provided by the host layer (substrates). Because of the complex nature of
this system exact solutions are not possible and we have used Molecular Dynamics

(MD) simulations to probe their physical properties (see below).

In chapter 4, the structure of domains and domain walls in the low temperature
solid phase of stage-n(2 2) alkali-metal graphite intercalation compounds (GICs) is
studied. For these systems, one can, to a good first approximation, ignore the effects
of interlayer intercalant interactions. Since detailed X-ray diffraction measurements
have been done in K, Rb, and Cs stage-2 G103 and excellent potential models are
available for RbCu, we focus on this particular system. Our aim is to see whether the
low temperature structures obtained from MD simulations can provide a consistent
explanation of the in-plane X-ray diffraction measurements. Based on the simulation
results we also want to see if simple Periodic Domain Wall (PDW) models can be

constructed to explain the dominant features of the X-ray diffraction experiments

in these GICs.

Chapter 5 is on the melting and freezing transitions of a 2D system consisting of
charged particles interacting with a repulsive screened Coulomb (Yukawa) potential.
The motivation for this work is two-fold. One is to compare our MD results with the
experimental data in R5024 and the other to study in general the role of corrugation
on the melting transitions of a 2D solid. In particular, we have investigated the role
of an incommensurate substrate corrugation potential of six-fold symmetry on the
melting transitions by probing the temperature dependence of the bond orientational
order parameter and its conjugate susceptibility. Other physical quantities such as
energy, diffusion constant, and the density of local topological defects have been
monitored through the transition region to find out the detailed nature of the phase

transition.

In chapter 6, we discuss the dynamics in the liquid phase of a 2D screened
Coulomb system by varying the strength of the corrugation potential. Mean square
displacement, velocity-auto correlation function, intermediate scattering function,
and dynamic structure factor are examined. Neutron scattering studies of the dy-
namic structure factor at large wave numbers from a 2D fluid on a corrugated sub-
strate show an extremely narrow resolution limited central peak on top of a broad
central peak, and a broad finite frequency peak. We have attempted to reproduce all
these features using MD simulations. Furthermore, since we can tune the strength
of the corrugation potential, we study the evolution of these spectral features as a

function of the corrugation strength and scattering wave number.

Finally, in chapter 7, we discuss the effects of corrugation on the phonon density
of states (PDOS) associated with in-plane vibrational excitations of the intercalants
in the solid phase. At very low temperatures, in the absence of corrugation, the pri-

mary excitations are phonons of a triangular lattice. In the presence of corrugation,

the ground state structure of Rng., can be approximated by a PDW solid with a
few “local” defects. The excitations of pure PDW solid which has a commensurate
ground state structure have a gap determined by the strength of the corrugation
potential. One important question we want to address is how the PDOS changes
when we go from the commensurate to the incommensurate system. PDOS associ-
ated with the in-plane vibrational dynamics of intercalants in corrugated 2D systems
such as GICs obtained from inelastic neutron scattering (INS) measurements at large
wave numbers show a characteristic two-peak structure. This two—peak structure
was explained previously by one phonon density of states, the phonon frequencies
having been obtained from an unscreened Coulomb interaction model. Earlier
MD simulation studies also supported this picture but with a screened Coulomb
model. The fact that the two completely different interaction models could give
the same PDOS implies that something is basically wrong with these explanations.
Using MD simulations and screened Coulomb model which has been extremely
successful for alkali metal GICs, we propose to clear up this confusion and provide

a correct understanding of the INS experiments.

Chapter 2

Structure of Disorderd Interfaces
in Coherent Multilayers

2.1 Introduction

With molecular beam epitaxy (MBE), one can routinely grow coherent strained su-
perlattices (CSLs) with I layers of A atoms and m layers of B atoms repeated 1) times,
denoted by (AiBm ), [1]. In thin (small I and m) epitaxial layers the lattice mismatch
between pure A and B layers is commonly accommodated by distortion of the unit
cells, resulting in strained pseudomorphic structures [2, 3]. Some of the superlat-
tices of current interest are lattice mismatched heterostructures for electronic and
optoelectronic device fabrications, for example (SizGem )p [4] and I nP/ I n1_,Ga,As
[5] superlattices. In the latter system, one can in fact control the lattice mismatch
and hence the coherency strain by controlling :c, a: = 0.47 giving the perfect lattice

match with I nP.

For ideal devices one should not have interlayer diffusion across the interface (to
avoid interfacial disorder) and misfit dislocations [6]. The misfit dislocations can be
avoided by growing thin superlattices but it is generally difficult to avoid interfacial

disorder [7]. Recent direct imaging studies of interfacial structure in an ultrathin

10

(Si4Geg)24 superlattice based on Z-contrast transmission electron microscopy [4, 8]
show that Si-layers contain a large concentration of Ge atoms. In I nP/ I n1-,Ga,As
[5] superlattices it is in fact possible to introduce interlayer diffusion (hence inter-
facial disorder) of In and Ga ions in the presence of Zn ions. Interfacial disorder
(IFD) in lattice mismatched systems leads to bond length fluctuations and hence
destroys the ideal physical characteristics of the device. Also under growth con-
ditions if IFD reduces the strain, then a disordered interface will be energetically
favorable. Thus it is important to understand the effects of IFD on the structural
properties of superlattice systems. In addition to the semiconducting superlattices,
there has been also a great deal of interest in recent years in magnetic multilayers
showing giant magnetoresistance (GMR) (Pratt at al. [9] and Levy et a1. [10]). It is
believed that the observed GMR [10] is sensitive to IFD. In this chapter, we address
the questions which deal with the structure of disordered interfaces by studying the
effects of IF D in (AgBm)p superlattices [11] on various physical characteristics such
as (1) the average bond lengths both parallel and perpendicular to the growth di-
rection, (2) specific bond lengths such as AA, AB, and BB both near the interfaces
and away from them, and (3) change in the coherency strain energy with degree of

IFD.

2.2 The Model

The model we have studied is illustrated in Figure 2.1. The bottom panel is a Z-
contrast scanning transmission electron microscopy image of a nominal (Si4Geg)24
superlattice [4] grown along [110] direction showing 32' and Ge multilayers. The mid-
dle panel is a schematic picture of CSLs with disordered interfaces. The white and

solid circles stand for two different types of atomic layers without mixing and the

11

pm /
7

 
   
  
  

 
 
   

       

  
   
         
       
     
  
 

' ye p. .u q s '
'. ...".".:H W......:.

       
 

  

C C-
C C C C. C C C
--C C .C C C C .C C'
C C C C'. C
--: .‘C C C C. :C C'
C C “C C
.. ... ..‘... .:. .‘... .

     

       

lei ii. I .

Figure 2.1: Heterostructures with interfacial disorder: The bottom panel is Z-
contrast scanning transmission electron microscopy image of a nominal (Si4Ge8)g4
superlattice [4]. The middle panel is a schematic picture of disordered interfaces, the
shaded circles representing the disordered layers. The top panel shows a continuous

profile p(z) of the concentration of one type of atom; w is a measure of the width
of the disordered region.

12

shaded circles represent the disordered layers whose layer stoichiometry depends on
the degree of mixing of two types of atoms. In the top panel we plot the concentra-
tion profile of one type of atoms, 1:”, for different layers (,1) treated as a continuous

function p(a:); w is a measure of the width of the disordered region.

For simplicity we will consider a 2-dimensional superlattice although our results
are general and can be applied to the 3-dimensional case. The energy associated
with the distortion of the individual bonds from their natural lengths L2,- is given
by

v = $2 KijULijl - L23)”, (2.1)
where Kg,- can be KAMKBB, and K43 and L?)- can be LangB, and Lag [12]. The
summation goes over all nearest neighbor bonds ij and the angular brackets imply
that the nearest neighbor sites are only counted once. We introduce a reference
lattice (I...) that is a perfectly undistorted triangular network and assume a small

displacement 11; associated with the site i.
Lu“ = I: + [(ua - u,-) ' EAR... (2-2)
where f. can take three possible values depending on i j and Ej’s are unit vectors

along the undistorted triangles of the reference lattice. By using harmonic approxi-

mation and minimizing the energy (Eq. (2.1)), we find that

z: K;,-[(u.- - Uj) ' EAR: - X1: Kij(L?jRa - 15) = 0- (23)
Eq. (21;) can be rewritten as follows:

A? -u = F, (2.4)

where F = K,j(L?jR,j — 1:). A? is a random matrix and F is a column vector with
random elements describing the deviation of the natural bond lengths from that of

the reference lattice.

13

When K;,- = K (i.e. when all the nearest neighbor interactions are the same),

It? becomes a nonrandom matrix and
u =51: (2.5)

where a is the usual Green function for the perfect system. In Eq. (2.5) disorder
appears only through the natural bond lengths L91- occurring in F. Knowing u we
can obtain the bond lengths ng,

Lij = i + K EXEC». - 5,2.) ' (RnkLgk - il- (2-6)
km

Now, we can calculate average bond lengths, their fluctuation, and strain energy.

2.3 Analytic Solutions and Numerical Simula-
tions

Up to this point, all results are formal and can be applied to any system. Now, we
will describe how to interpret Eq. (2.6) for the CSLs with a specific lattice structure
i.e. triangular lattice. First, I: has two distinct directions; one (6) is along the layers
and the other (V or n) connects neighboring layers. Second, each layer has different
average concentrations of A and B atoms. Here, we use Greek letters to denote the
layer index. Therefore, '1' can be rewritten as pn which corresponds to the nth atom
in the ,uth layer. Concentration of large atoms B in the layer p is denoted by 1:“.
The layer concentration of small atoms A is therefore 1 — m”. For the CSLs with

ideal interfaces, 2:“ = 1(0) for the layers with all B(A) atoms.

In the absence of IFD and in the limit of uniform force constants, i.e. K .3 = K,
I: can be determined by demanding that there is no net macroscopic force on an

imaginary line through the reference lattice. From Eq. (2.3) when K ,3 = K,

14
K:(L— L? Rs)- - (2-7)

Therefore the nearest neighbor distance of the underlying virtual crystal for an

arbitrary direction is obtained from the above equation.

We introduce a projection operator 0",, to deal with the random variable L2,.
0“,, = 1(—1) for an A(B) atom in the pth layer. The mean value of 0,". for a given

layer p is
()0” =Z0un/N- — l — 22:”, (2.8)

where N is the number of sites in the layer. For the sake of simplicity we will use
only the layer index p instead of pn to describe 0,". since (0“,) is independent of
the site index n. However, when the random variable 0,, appears before taking an
average, its 12 dependence is implicitly assumed. Consider a nearest neighbor bond

between atoms in two layers p and 6. The natural bond length is given by

 

 

 

 

I 1+01+0€ 1—01-0:
2'5 2 u 2 9” 2 “ 03
1+ 1— 1— 1+
+ [ 0,, 20£+ 20" 203] L03 (2.9)

If we take A E L?“ + L933 — 2LAB = 0 i.e. L943 is the arithmetic mean of L9“ and

L993, then Eq. (2.9) reduces to

1+0 1+0 1—0 1—0
L3“: ( 4“+ 4‘)L3A+( 4“+ 4‘)Lg3. (2.10)

An important difference between Eq. (2.9) and Eq. (2.10) is that they lead to differ-

 

 

ent reference lattices; the former gives an isosceles triangle whereas the latter gives
an equilateral one. As a result, when A = 0, we need only one topological parameter
to describe the rigidity of the reference system which is an equilateral triangle. Also
the resulting equations for different physical quantities are much simpler. From here
on we will discuss our results only for the case A = 0 and results for A 7E 0 will be

given in Appendix A.

15

2.3.1 Average Bond Lengths

First, we calculate the average distance (14,, ,g) for inter- and intra—layer spacing
denoted as V (or n) and 6 directions, respectively. We find for the average intralayer

bond length
(Lu.u>5 = (1' X)L?AA + XLgBa (2'11)

where X is the global concentration of B atoms over the entire CSLs i.e. X =

2;, 3,, /N1, where N; is the total number of layers. Along u (or 1]) direction
”flatly = (1 — $u.€)LOAA + xu.£L%Ba (2°12)

where

xuvfi = 3,1341%. (2.13)

 

Thus the average bond length along the layer (coherency strain direction) depends
only on the global concentration X of B atoms due to the requirement of coherency
whereas that between the layers depend on the local concentration. This asymmetry
is a characteristic of the CSL structure and is not found in random alloy systems

[12].

In the absence of IFD, :12“ has a steplike structure, 0 in A multilayers and l in B
multilayers whereas in the presence of IF D this steplike behavior of 33,, changes to a
smooth function of )4 near the interfaces (see the top panel of Figure 2.1). Numerical
simulations with conjugate gradient method have been done for the system (.4le )p
with I = m = p = 5 and with 200 atoms in each layer (see Figure 2.2 where one unit
of A585 is shown). We have chosen L3“ = 1 and L933 = 1.04; a 4% lattice mismatch
is typical for semiconductor alloys. Each layer can be denoted as .41-,"3," with

3:“ changing from one layer to the other. For simulation, we considered IFD to be

16

 

 

 

 

Figure 2.2: A schematic picture of a coherent superlattice with disordered interfaces.
The multilayer structure is (AsBs)5. Layers 2, 3, and 4 have all A atoms. Layers 7,
8, and 9 have all B. atoms. Layers 5, 6 and 10, 1 are interface layers, shaded circles
representing the disordered layers. 6, V, and 1] are directional vectors we use for the
reference lattice.

17

confined to the two interfacial layers only. In this case, x“ = 0 for p = 2,3,4 and
z, = 1 for p = 7,8,9. The 5th and the 6th layers (as well as the lst and the
10th layers) are interface layers and we have considered a situation where the A
atoms from the 5th layer and B atoms from the 6th layer are mixed but the total
number of A and B atoms are conserved within the two adjacent layers. Because of
the conservation of the number of each type of atoms, we have 1 — $5 + l — 1:6 = 1
(2:54-36 = l) for A (B) atoms respectively. Thus the 5th (6th) layer forms a random
alloy 141-:5st (Ar-2.3:. E ArgBl-:5)' I

Figure 2.3 gives the projection of average interlayer spacings (Lumfi), along the
layer growth direction (denoted as y) as a function of the layer index it for several
profiles of IFD i.e. from no interlayer diffusion (perfectly clean interface, 2:5 = 0 and
23 = 1) to complete exchange (layer interchange, 2:5 = l and $3 = 0). The inside
four different curves are for different degrees of interface mixing i.e. 3:5 = 0.1, 0.3,
0.5, and 0.7. In all cases, numerical simulations (denoted by symbols) and analytic

calculations (solid lines) agree very well.

To see how the average interlayer bond length (or equivalently the average bond
length joining atoms in the layers [1 and p + 1) changes with IFD, let us look at
Figure 2.4. The lowest ((L3,4)) and the tOpmost ((L7,8)) lines are independent of
2:5 simply due to the absence of any interlayer mixing. Similarly, (L253) and (L83)
are also independent of $5 and have not been shown in the figure. The linear
increase (decrease) of (LN) ( (La-1)) can be easily understood from Eqs. (2.12) and
(2.13) by noting that $4 = 0,1:-; = l, and 2:6 = l — 2:5. One interesting result is
that the average bond length between layers 5 and 6 (L55) does not depend on
the concentration 1:5 although these are interface layers and the concentration of B

atoms in each layer changes. This is because the average bond length (Eq. (2.12))

18

 

  

 

 

 

 

0.91 L—I I I T I I I I I I T UWT—l V I l I I U V U I_.:
0.90 '— —;
r .
r- u
h0.09 -_—-
’E.
*i .
0.88 —
A . ..
v x - .
0 ; -
0.87 D 3 -
+ : Z
X : 3
0.86 _ a : x5 = 1 ‘.
l l l l J 1 l 1 l I l l l l [J I I l l 1 L1
0 2 4 8 8 10

,u.

Figure 2.3: The same system as Figure 2.2. Comparison between analytic calcu-
lations (solid lines) and numerical simulations (symbols) of average layer spacing
(mel), as a function of the layer index p for the multilayer (A585)5. 15:0 means
perfect interface whereas 2:5 = 1 implies the interchange of two adjacent interface
layers.

19

 

 

 

 

 

 

 

1.05 -T I I I I I I I ITIT I I l I I I I I I I I I d
I 1

1.04 x x x a: x
1.03 - _
’3. r -
51 02 ' v v .. v v '
.3 I 3
V _ .
1 01 — —

1.00 X x x x x
C 2

0.99 l l l l l l J l l l l l l I 1 l l l I l l l 1
00 02 0.4 06 08 1 0
x5

Figure 2.4: The same system as Figure 2.2. Each line corresponds to an average
bond length (LAHI) between layer p and p + 1, 2:5 is the concentration of B atoms
in the layer 5. From the bottom to top p is 3, 4, 5, 6, and 7. The x symbols are
numerical simulation results.

20

depends on the average concentration (Eq. (2.13)) of the two layers and since the
change in concentration of the atoms occurs only through the exchange of A and B
atoms within the two layers 235,6 = .25 + $3 = 1. Therefore, there is no change in
the average concentration (2:55) between the interface layers resulting in (L55) being
independent of IFD. Each layer follows Vegard’s law [13] and numerical simulations

and analytic calculations agree very well.

2.3.2 Specific Bond Lengths

We have calculated specific bond lengths both along (denoted as 6-bond) and be-

tween (1/ bond) the layers. For the 6-bonds (same layer), we find

(Lust 3“ = (E) + a":c,,(LgA " L238) (2'14)
(Lumlsaa = (I’) —' an“ " $u)(L?4A — [1238) (2'15)
(L... 3‘” = (L) — aw; — 2:009... —- Lass). (2.16)

Here, the mean length (L) = (l — X )Lg A + X L93 3 is determined only by the global
concentration X due to the coherency requirement along the 6-direction. The pa-
rameter a" [14, 15, 16] in the above equations characterizing the topological rigidity

of the reference lattice is given by
It 1 - fl - I - I
a = gZESIMQ-flh- D (Q) l"rlsnrl(Q-7), (2-17)
71’ Q

where z is the number of the nearest neighbor atoms and 7 and 7’ stand for the
nearest neighbor directional vectors (16, ill and in) of the triangle describing the
reference lattice. It is easy to see that (Z) = (1 — :1:,,)2(L,,,,,)g‘A + xfi<Lu,u)5BB +
2:0,,(1 — zu)(L“,u)§B. Furthermore, because of the assumption A = 0, we have an

additional relation (LN),M + (Lu,,,)6BB — 2(Lum)?‘9 = 0,

21

Figure 2.5 gives the 3,, dependence of different normalized intralayer spacings
for the layers 5 and 6. When 23;, = 0, (L,,-;,5);5‘A = 06,6)?” and these depend only
on (L). Our choice of I = m = 5 gives the global concentration X = 0.5. Therefore
(L) is 0.5 in the normalized length. Here, we introduce a normalized bond length
such that

L AA __ L0
(du£>AA = ( £§>5 L0 AA. (218)
BB — AA

When :05 = 0, the AA(BB) bonds in the 5th (6th) layer are already expanded
(contracted) to form a CSL. Also, (d“.,,)3“((du,u)5BB) is never equal to 0(1) which
reflects a characteristic (or requirement) of the CSL as we see in Figure 2.5. Let us
try to understand the $5 dependence of ((155); in Figure 2.5. When B atoms from
the 6th layer go into the 5th layer making $5 96 0, the stretched AA bonds tend
to go towards their natural bond lengths L9”. This is why ((15.5); decreases as $5
increases. On the other hand, when $5 = 0, the layer 6 (A1_,,B,, E A,,B1_,5) is
filled with all B atoms which are contracted to match with layers with A atoms. As
:05 increases i.e. atoms with small size come into this layer, large atoms can relax

toward their natural lengths, causing an increase in (d6,6)5 with 1:5.

For the u—bonds (interlayer bonds) connecting two adjacent layers with concen-

trations 3,, and 2:5, we obtain

(Muslim = (Lac). +a"$u.L°£( AA - L0 as) (2.19)

(Luxlfa = (Lac). ‘ 0"(1 - mu.<)(L3A - Lies), (2.20)
.. 1

(Ll-‘95):B = <L“v£>u — a (5 — $Hv£)(L91A — L%B)’ (2'21)

where (10%)., and 27“,: have been defined in Eqs. (2.12) and (2.13), respectively. The
average interlayer bond lengths are also determined by the same single topological

rigidity parameter a”. In the V or 7] direction, when 104,5 = 0 i.e. when :05 = 0, the

22

 

1.0IIIIIIIIII'IIIIIIIIIIIIfI

0.8

0.8

<d6.6>6

0.4
0.2
0.0
0.8
0.8

0.4

<d5.5>6

0.2

 

0.0lljlllllillJlllllllllllj
0.0 0.2 0.4 0.6 0.8 1.0

x5

Figure 2.5: Normalized average intralayer bond length for the interface layers as a
function of :05. Three curves of each set from top to bottom correspond to BB(<>),

AB(><), and AA(D) bond lengths.

23

AA bond is unstrained (see Figure 2.6). Here, 24 is always equal to 0 because we are
considering the interface mixing only between the layers 5 and 6, i.e. the layer 4 has
all A atoms. Therefore, possible bonds between the layers 4 and 5 are only of the
AA and AB type. Similarly, the layer 7 has all B atoms and bonds between layers
6 and 7 are only of BB and AB type. As B(A) atoms come into the layer with all
A(B) atoms, interlayer bond length ((0.5),. ((d6,7),,) becomes bigger (smaller) which
means that the AA(BB) bond expands (contracts). Therefore, “45);?“ increases
as the concentration of B atoms in the 5th layer increases. The $4.5 dependence of
average AB bond length can be explained similarly. Of course, all combinations are
possible across the layers 5 and 6. However, ((155),, for AA, AB, and BB bonds are
independent of the 35,6 and can be explained by the same reasoning as discussed for

the Figure 2.4 in the previous section.
2.3.3 Fluctuations and Strain Energy

Next, we calculate the global bond length fluctuations and strain energies associated
with the bonds in the 6- and V-directions, respectively. In particular, we study how
the strain energy depends on the degree of IFD and a" and investigate under what

conditions strain energy favors IFD.

For fluctuations of bond length distribution, we find

2 _ 2: u _ m- o _ 0 2l+mxu(1—xu)
(L )6 (L>6 a (1 a )(LAA LEE) 2: 2(l+m) , (2'22)

(L2). - (L)3=<L2)n-(L>2

7)
l+m

0"(1 - a")(L?u - L933)? 2
u—l

$141— 3:4) + mil-10 _ mil—1)
4(1 + m)

 

. (2.23)

We see that bond length fluctuations over the entire CSLs are same for the 6 and
:2 directions in spite of the coherency constraint along the 6-direction. The specific

bond length fluctuations for AA, BB, and AB are also found to be the same. The

24

 

 

 

 

 

 

 

 

 

 

 

1.0:
0.8:
A: :
1.. 0.6:—
0 I
“U 0.4.— '3.
v 5 :
0.2:- 1
:lllllllllllllllllllllJll:
°,0E IYIIIIIIIIIIIIfiIIIIUIrIE
)- III
0.8:- '1
A; §——¢ 0 o——§
o 0.8:- a
,5 a x x s:
: c4
'6 04.— a fi-:
V a U U -I
0.25- —:
00::ili:f%§§:l§§§:ili§:lift:
E E
0.8..— '1
A: : i
,q 06:— '1
* :
“U 0.4:-
V E
0.2:
: lllLll IlLl
0.0 1
0.0 0.2 0.4 0.6 0.8 1.0

Figure 2.6: Specific bond lengths near the interfaces for u-direction vs. :05 concen-

tration of B atoms in the layer 5. The symbols stand for the same as Figure 2.5 i.e.
BB(o), AB(x), and AA(D) bond lengths.

25

average strain energy per bond is given by

 

 

0 _ 0 2

5 = W“); + 5., + 8"), (2.24)
where

1 l+m 2 to z#(1 — 3”)
85—ng (X—zu)+(1—a ) 2 , (2.25)
__ In (+171 _ _
a, = a, = Q—i-l “(I 3,) + “"0 0‘"). (2.26)
1+ m “=1 4
Combining Eqs. (2.24)-(2.26), we obtain
1 (+77: 13

6 = €09... — L25)“ $0 - a)? + JET-1&0 — an] . (2.27)

I + m ”:1
where 2c is the coordination number in the 6 direction (along the layer) and z is the

number of nearest neighbors.

If we describe the interfacial diffusion by a continuous function p(:z:) of width w

(see Figure 2.1), we obtain

= 52:02... — L13)” 3m”) - x2) + 5;;90 — am] , (2.28)

where (3") = f z"p(z)d:r. Thus the difference in the strain energy between the
multilayer with perfect interface (w = 0) and that with IFD is given by

261

a.-. — a. = 4210?... — L28)” ; — -,-(l — a")] (($0.20 — 02).). (2.29)

Since the term (x2)w=o — (3:2),” is always positive, we see that when 3:- > %(1 — a"),
Ew=o > 8,” i.e. a perfect interface has higher energy than the disordered interface.
Therefore, in this case, interface mixing is preferred due to the reduction in the
strain energy. For a triangular lattice, a" = 0.392, 2 = 6, and 2c = 2. Therefore the
IFD will reduce the strain energy of this system. For the FCC lattice, a“ = 0.24

and z = 12. In this case if the growth direction is (100) then 2, = 4 and a disordered

26

interface increases the strain energy thereby forming a clean interface. However, if
the sample is grown along (111) direction, then 26 = 6 and in this case interface
mixing will be preferred. This can be understood easily because 26 measures the
number of bonds in the coherency plane. The larger value of 2,; results in more

strain energy because more bonds are constrained to make CSLs.

Before we close this section, let us discuss the result for the strain energy when
A 9% 0. For the sake of simplicity, we have again considered interfacial mixing of the
two adjacent layers 5 and 6 only, the same situation for which numerical simulations

were done. Total strain energy is
Eta! = 5 + 5A, (2.30)

where 8;; is the additional term contributing to the strain energy due to nonzero A

which is as follows.

K

8A = 325(1 — $5)A2 [3 ._ 63:5 + 22:: _ m

1+ m
— 2(1 " “X13035 + (2 - Islafm) - 2(1 — 215)((1 — 235W; + 25:0)

— 4x5(l — 22:5)(03, + 53,) -4x§bf,,,] , (2.31)

where the parameters 036, a;,,, 036 etc. are defined in the Appendix A. Figure 2.7
shows the dramatic difference in the x5 dependence of the strain energy for various
values of A for the triangular net. For A = 0, the strain energy vs. concentration
of the interface layer shows a minimum at .25 = 0.5 (the bottom curve), consistent

with our prediction. However, as we increase A, IFD is suppressed.

2.4 Summary and Conclusion

We have investigated the effects of interfacial disorder on the structure of coherent

multilayers of two types of atoms. A simple central-force model has been set up for

27

 

0.15IIfiIIIIIIIIIIIIIIIjIfiII

0.14

Em

0.13

 

0.12

D

 

 

llllllllllllllllllllllll

 

0.0 0.2 0.4 0.6 0.8 1.0
x5

Figure 2.7: The strain energy vs. concentration of the interface layer 5 for the CSL.
The system for the A = 0 (bottom curve) shows the minimum strain energy at
$5 = 0.5 which indicates that interface mixing is preferred. Here, A = 0 means that
L1, = 1, L33 B = 1.04, and LAB = 1.02. The other curves are obtained by gradually
changing L33 with L33 = 1.015, 1.01, 1.005, and 1, the top curve corresponding to
L213 = 1.

28

this purpose. To preserve the coherent structure, the size of the system parallel to
the layer is fixed. The problem of disordered interfaces can be solved exactly in the
harmonic limit and when the spring constants between AA, BB, and AB atoms
are the same. We find that the average interlayer spacings, the strain energy and
average bond lengths depend on a single topological parameter a“ when A = 0.
Layer spacings change in the presence of interfacial disorder and can in principle be
measured by X-ray experiments. Individual bond lengths have been calculated and
can be obtained from EXAFS and nuclear quadrupole resonance experiments which
are in practice difficult because of the small signals from the interface regions [17].
We also find that the strain energy can determine the degree of interfacial mixing
when L9“, L933, and L943 are given although other sources such as charge transfer

energy may be relevant in some situations.

Bibliography

[1] L. Esaki and R. Tsu, IBM J. Res. Develop. 14, 61 (1970).

[2] DB. McWhan, Synthetic Modulated Structures, Ed. by Leroy L. Chang and
BC. Giessen, Academic Press, New York (1985), p43.

[3] An example of a strained layer superlattice in InGaAs/InP system has been dis-
cussed by J .M. Vanberg, D. Gershoni, R.A. Hamm, M.B. Parrish and H. Temkin,
J. Appl. Phys. 66, 4035 (1989).

[4] DE. Jesson, SJ. Pennycook, and J. M. Baribeau, Phys. Rev. Lett. 66, 750
(1991).

[5] D.M. Hwang, S.A. Schwarz, T.S. Ravi, R. Bhat, and CY. Chen, Phys. Rev.
Lett., 66, 739 (1991).

[6] J.W. Matthews and A.E. Blakeslee, J. Cryst. Growth 27, 118 (1974); ibid 29,
273 (1975).

[7] A. Bourret, P. Fuoss, G. Feuillet, and S. Tatarenko, Phys. Rev. Lett., 70, 311
(1993).

[8] S.J. Pennycook and DE. Jesson, Phys. Rev. Lett. 64, 938 (1990).

[9] W.P. Pratt, Jr., S.-F. Lee, J .M. Slaughter, R. Loloee, P.A. Schroder, and .1. Bass,

Phys. Rev. Lett. 66,3060 (1991 ).
29

30

[10] Horacio E. Camblong and Peter M. Levy, Phys. Rev. Lett. 89 2835 (1992).

[11] Hyangsuk Seong, S.D. Mahanti, and M.F. Thorpe, Bull. Am. Phys. Soc. 36,
991 (1991); preprint.

[12] M.F. Thorpe and E.J. Garboczi, Phys. Rev. B 42, 8805 (1990).

[13] L. Vegard, Z. Phys. 5, 17 (1921).

[14] Y. Cai and M.F. Thorpe, Phys. Rev. B 46, 15872 (1992).

[15] Y. Cai and M.F. Thorpe, Phys. Rev. B 46, 15879 (1992).

[16] Normand Mousseau and M.F. Thorpe, Phys. Rev. B 46, 15887 (1992).

[17] Z.H. Ming, A. Krol, Y.L. Soo, Y.H. Kao, Bull. Am. Phys. Soc. 38, 84 (1993).

Chapter 3

Gallery Expansion of Randomly
Intercalated Anharmonic Bilayers

3.1 Introduction

All crystalline solid solutions show a composition dependence of the lattice constant
which increases with the concentration of the large component (3:) [l]. The well-
known Vegard’s law [2], where the average volume depends linearly on x, is obeyed a
large class of semiconducting and insulating alloys of the type of A1_,,B,C [3]. How-
ever, it is not obeyed by metallic alloys and ternary intercalation systems. These
systems exhibit rather complex nonlinear behaviors which depend on the compe-
tition between local and global energies associated with the formation of a solid
solution. An important class of such solid solutions with which we are concerned
in this thesis are ternary intercalated layered systems, A1_,B,L with 0 S a: S 1,
where A (small) and B (large) are two different intercalants and L represents the
host layer. The major expansion of these materials occurs along the direction normal

to the layers, which is commonly denoted as the c-axis[4]-[9].

Several theoretical models [5],[10]-[l5] have been proposed to understand the

nonlinear z-dependence of the c-axis expansion and the condition under which Ve-

31

32

gard’s law can be seen. In one limiting version of these models which has been
referred to as the “rigid-layer model” [5, 10], the host layers are assumed to be flat,
i.e., completely rigid against transverse distortions. In this limit, the nonlinearity of
the c-axis expansion arises from the finite but different compressibilities of the two
intercalants i.e. different strengths of the harmonic interaction between the inter-
calants and the host. In another limiting version [10]-[13], so called “layer rigidity
model”, the finite transverse layer rigidity of the host layers has been taken into ac-
count but the intercalants have been treated to be completely incompressible (hard
spheres of different diameters) and any finite concentration of large spheres opens up
the gallery to its maximum value. In this model, the nonlinear z-dependence arises
from the transverse rigidity of the host layers. To circumvent these two extreme as-
sumptions, namely infinite layer rigidity and completely incompressible intercalants,
a simple harmonic spring model [14, 15] has been introduced. In this model, the
host-intercalant interaction is represented by two harmonic springs of equilibrium
lengths ha and 11% with spring constants K A and K B, respectively and the host
layer rigidity is characterized by a spring constant KT. This model has been in-
vestigated thoroughly both for bilayers [14] and for multilayers with and without
correlation between the occupation of the intercalants [15]. However, this harmonic
model, although quite general in many respects, does not address the role of anhar-
monic interaction in the gallery expansion problem. In some cases, it does not give
a proper account of experimental data near a: = 1 [11] and the discrepancy becomes

large for rigid layers and nearly incompressible intercalants.

On physical grounds one can expect the anharmonicity in host-intercalant in-
teraction to be important when the layers possess finite rigidity against transverse

distortion. In Figure 3.1, we show a small intercalant surrounded by large inter-

F_'P

33

calants. If the host layer is completely floppy then one essentially samples the
harmonic region of the host-intercalant interaction potentials both for the small
and large intercalants. On the other hand, if the host layer is rigid, although the
interaction between the host layer and the large intercalants can be treated in har-
monic approximation, one samples the anharmonic region of the potential between
the small intercalant and the host layer. Harmonic approximation is not appropriate

to model this latter interaction.

In addition to the atomic potential models described above, a discrete host layer
distortion model with hard sphere intercalants has been proposed [10] and solved
exactly [12], and its results have been applied to fit the experimental data [11, 13]. In
this model referred to as “Catchment Area model” (CAM), the average normalized
interlayer spacing is given by (d) = 1 — (1 — 2:)", where a: is the concentration of
the large intercalant and p is a measure of the healing length. In other words, p
is a measure of the size of the catchment area i.e. the region over which a large
intercalant B with size hg opens up the gallery spacing in its neighborhood to hg.
Also, in this model, it is assumed that if another large intercalant is present within
this “healing length”, then the gallery spacing at this site does not increase any
further and remains at h%. The form of (d) obtained with CAM describes the gallery
expansion of layered solids such as “1.1/1-sz (.V is a vacancy) and Cslethm
(Vm is a vermiculite) [13] very well. From the fit to the experimental data, one finds
that the host layer rigidity parameter (p) increases from graphite to vermiculite and

is consistent with the increased transverse rigidity of the host layers.

The CAM as described above represents an extremely anharmonic system, as
regards both the host-intercalant interaction and the host layer distortion are con-

cerned. Unfortunately, in this model anharmonicity cannot be controlled. It is

"i

34

(a)

@002
o

9 013)

(C)

 

 

Figure 3.1: An example of a 1-dimensional bilayer model in which a small intercalant
is surrounded by large intercalants. The solid lines represent the host layers and
the circles depict the intercalants. The diameter of the circle A(B) is 1130123) which
is the natural length of atom A(B). (a) The host layer is perfectly rigid; therefore,
equilibrium length of A atom, hA, is h%, thus sampling the anharmonic region of
the potential between the small intercalant and the host layers. (b) The host layer
is perfectly floppy; therefore, h A is ha which is its natural length. (c) The host layer
has a finite rigidity; therefore, ha < hA < kg.

35

therefore important to study a physical model where one can control the degree
of anharmonicity and see how it affects the average gallery spacing and the local

gallery distortion.

In this chapter, we introduce an anharmonic potential model and study this
using both numerical simulations and analytic techniques. For simplicity, we use a
bilayer model, in which we only consider a single gallery between two host layers
where the intercalants go. Also, we include anharmonicity only in the potential
energy between the host layer and the intercalants, while using the usual harmonic
approximation for the host layer deformation energy. Preliminary results of this
model have been reported in Ref.[l6]. Taking account of anharmonicity effects in
the host layer distortion energy is extremely difficult and will be the subject of future

research.

The arrangement of this chapter is as follows. In Sec. 3.2, different anharmonic
potentials between the host layer and the intercalants are introduced and the results
of a few limiting cases (perfectly floppy and perfectly rigid layer) are discussed.
In Sec. 3.3, the single impurity problem is studied analytically to understand the
behavior of the system near a: = 1. An effective medium theory is developed in
Sec. 3.4 and its results compared with numerical simulations in Sec. 3.5. In Sec. 3.6,
we compare the results of the numerical simulations of Lennard-Jones potential with
that of the CAM discussed earlier. Finally, in Sec. 3.7, we give a brief summary of

our major findings.

3.2 The Model and Limiting Cases

Consider a ternary layered system A1_,B,L consisting of a single gallery with inter-

calants randomly occupying a set of lattice sites between two host layers. A (small)

36

and B (large) are two different types of intercalants and L is the host layer such as
graphite, dichalcogenide, or vermiculite. For a fixed intercalant distribution, there
are two major contributions to the total deformation energy of the system; one is
the interaction between the host and the intercalants and the other between the

host atoms themselves,
1
E = E11] + '2- ZKTUln — hn+l)2- (3.1)

In Eq. (3.1), Em is the interaction energy between the host and the intercalants,
K1- represents the transverse rigidity of the layers, and h“ is the local gallery height
(the distance between two adjacent host layers) at the n-th site. In the harmonic

approximation, E”; is replaced by harmonic potentials [14], i.e.,
1
EH! = 5 z Kn(hn - ’13.)“, (3.2)

where K, and h: are the spring constant and the equilibrium height of the inter-

calant at the site n, respectively.

A simple and yet realistic way to extend the model to physical systems is to
consider a Lennard-Jones (LJ) potential to represent the interaction between the
host and the intercalants. For charged intercalants, as in ternary complex layered
oxides, the interaction between the host and the intercalant is long ranged. But
for short wavelength “local” deformations, this long range interaction may not be
very important. In graphite intercalation compounds, the long range interaction is

usually screened and a short range model is adequate. In the LJ potential model,

an 12 0n s
Em = 31.1 = 246'. (—) - (—) , ' (3.3)
a hit hn
where the parameters 6,. and 0,, can be obtained, for example, from the gas-phase

data [17, 18]. To see the strength of the anharmonic effect in the LJ model, one

“a

37

can make use of Taylor series expansion of the LJ potential about the minimum

(h: = Wan)! namely,

 

 

 

1 1

EL, = Z": [13,, + -2—!K,,(h,, — 1:9,)2 — 57,01, — h:)3 + - - ] , (3.4)

where
6’13“ 726,.

K" ‘ aha, " {/203’ (3'5)

and
33E“ 6n
7,, 0h: — ”Rx/20: (3.6)

The second term on the right hand side of Eq. (3.4) corresponds to a harmonic
potential and the subsequent terms to the anharmonic contributions. In the LJ
potential, the strength of the first anharmonic term relative to the harmonic term
7,, / K n is fixed for a given value of 0... Therefore, to isolate the effect of anharmonic-
ity, we introduce a slightly modified model in which one can vary the strength of

the anharmonic term independently, namely,
1 1
EC 5 ‘27 g-Knaln _ hay - '3'? 227110“: _ hills, (37)

where Kn(> 0) and 7,,(> 0) serve as control parameters. In this model one can
easily study how the physical properties of the system evolve from the harmonic
interaction limit. This cubic potential model is simpler than LJ potential for analytic

calculations, hence we will concentrate on this model for analytic calculations.

In connection with the transverse rigidity of the host layer, one can think of
two limiting cases, KT —* 0 and KT -+ 00. If the layers are completely floppy to-
wards transverse distortions (K T = 0), then obviously all hn’s have their equilibrium

lengths. Therefore the average gallery height (h), defined by

 

38
(h) 2 +211... (3.3)

obeys Vegard’s law (i.e. (h) = (1 — z)hf’, + xhg) for both Eqs. (3.4) and (3.7). If
on the other hand, the layers are infinitely rigid (KT —+ 00), then the energy E in
Eq. (3.1) is minimized by having all hn’s equal. Then the average height has the
value

 

(1 — z)KAhfl“ + zxahgr (3.9)

h =
< l“ ((l-z)KAh‘},8+zK3hg8

for the LJ potential, and

(1 — z)(KA + 7.13) + z(Ka + 731:9.)
(1 - xl‘lx + 3‘78

 

(’00 =

 

_\m1— z)KA + zKB + 2:73AM“ — {(1 — z)“ + 273} Aha:(2K3 + 73Ah)
(1‘ 307.4 + $73 ’

 

(3.10)
for the cubic potential and Ah = h“); — 11%. In the harmonic limit (7,; = 73 = 0),
Eq. (3.10) recovers the results of Dahn et al. [5] which is

(1- $)KAhg + $K3h23
(1— $)KA + 23KB

 

(1).; = (3.11)

and obeys Vegard’s law when K A = K B.

We have performed a series of numerical simulations for both the LJ and the
cubic potential models. Typical behavior for the z-dependence of the normalized
spacing (d) is shown in Figure 3.2, where the simulation results for K A = K B are
compared with Vegard’s law. When K A = K 3, the harmonic potential system

shows Vegard’s law regardless of the value of K T [14]. However, even if K A = K 3,

 

39

      
   

 

 

1.
O '- I I l I l I I I FI I I T I I Ij I I J
3/
1- II
b d
0.8 - '-
L- ..
I- II
I- '4
0.6 - '-
+ 1
G . .
v , .
0.4 h— —
1- u
0.2 - '-
- Cl
, -I
l l 1 +1 I l l I] 1 l l l l l l 1 l 1 L1 L1

 

 

0.0
0.0 0.2 0.4 0.6 0.8 1.0
X

Figure 3.2: The deviation from the Vegard’s law for the anharmonic potentials:
K A/KB = 1,KT/K3 = 10. o : LJ potential, x : cubic potential with 7,, given by
Eq. (3.5), - : Vegard’s law. The solid lines on symbols are guide for eyes.

40

or equivalently 64/0} = 63/0123, (11)“ and (MC do not obey Vegard’s law but show

a nonlinear a: dependence due to the effect of anharmonicity.

The upward shift of the curves from the linear :c-dependence in Figure 3.2 can
be explained from the shape of the potentials used. For positive 7,,’s both forms of
the anharmonic potentials cost more (less) energy under compression (expansion)
than the harmonic potentials. Therefore, when one minimizes the total deformation
energy, the small atoms with the anharmonic potentials adjacent to big atoms tend
to stretch longer than those with the harmonic ones, while the big atoms tend
to compress less. This is manifested in the larger average gallery expansion for
the anharmonic potentials than that for the harmonic ones as seen in Figure 3.2.
If one had used negative 7,,’s, one would have seen downward shift (sublinear x-

dependence) of the curves from Vegard’s law behavior.

3.3 Single Impurity with Anharmonic Cubic Po-
tential

To understand the upward shift (or what is referred to as overhang) of the gallery
height quantitatively near a: = 1, let us consider a single small atom A in a chain of
N lattice sites with all other sites occupied by B atoms. For simplicity, we will use
a harmonic potential for the B atoms and a cubic anharmonic potential for the A
atom. Harmonic approximation for the B atoms has been found to be adequate from
numerical simulations. The anharmonicity of B atoms hardly affects the average
gallery expansion because large atoms are always compressed in bilayer systems.
Suppose that a single A atom is at the origin, then the deformation energy of the

chain is given by

41
Es=-23Kn..—(h 312,741... +5320.— hn+1)21 (3.12)

where

Kn = KB(1— 6m) + KAISno = KB - AK5no,

and
7n = 76m-

1.-.“!

Here, we have used AK 5 K 3 — K A and Ah 5 1193 — ha. Minimizing this energy

with respect to hn, we get
1 _
K,(h,, — hfi) — 57,01, - hi1)2 + KT(2h,, — h,+1 — h,_,) = 0. (3.13)

Next, taking the Fourier transform of h“ i.e. h" = 2:9 h(q)e"" and using periodic

boundary condition (hN+,, = ha), we obtain

(KB + 2KT(1 — cosp))h(p) — — [11:32 h(p- q) + —(K3Ah + 1° AK)

hO’ Ah , 1 , ,
7(- 3‘--hl’aZMp-qHTXZhW-q-qH522h(9)h(p-q-9)
9 9 q’ 0 q'
= KBh$6p’o. (3.14)

The quantity of physical interest, (h), which is equal to h(0) can be obtained from
the solution of Eq. (3.14). Assuming that 2,, h(p — q) = 2' h(q) and changing the

42
notation UN to 1 —- :0, we have the final result for the average gallery height,

(h) = h(°)=h%

       

W

 

 

 

1
+ 7W3.1_;€;[_AK_ \/{I:V-—B' - AI(}2 - 2K37Ah(-u7 - l) .
(3.15)
In the above equation, W is the Watson integral [14] given by
1 K” (3.16)

W=IV;KB+2KT(1 —cosp)

and is inversely proportional to the transverse layer rigidity KT. The results of
Eq. (3.15) for difl'erent K1 are compared with the results of numerical simulations
in Figure 3.3 and fit very well with the slopes of the average gallery spacings obtained

numerically.

When 7 = 0, the third term in Eq.(3.15) vanishes and the average gallery height
reduces to the correct harmonic limit of Ref. [14]. In addition, for 'y > 0, the
anharmonic contribution to (h) is always positive because 0 S W S 1. This positive
deviation results in a larger gallery expansion than the harmonic case caused by the
cubic anharmonic potential. Further, the more rigid the layer is, the larger this

deviation is.

3.4 Effective Medium Theory (EMT)

To understand the gallery expansion of A1_,,B,L for arbitrary values of 2:, we use
an effective medium theory (EMT) for the bilayer model with many defects. In this
section, the cubic potential is used for both A and B intercalants. To develop an
EMT for this anharmonic system, we follow the general ideas discussed in a series

of papers [14, 15, 19] developed for the harmonic potential model.

 

43

 

  

 

 

100 '- I I I l I I I I TM r I I q

.- D ..

_ a ..

0.8 — D ._

I 0 :

: U o :

0.6 — _

A : o X :
v I
V t x .
0.4 :- D 0 1

. x _

b X : Kr/Kn=o.01 -

0.2 — x 0 : K,/K,,=0.1 _

I D : K,/K,,=1.0 ;

o O._ l l l l l l l L l l l 1 I l l l l 1 I I l l I l q

 

'0.0 0.2 0.4 0.6 0.8
X

l"
0

Figure 3.3: Comparison of the numerical simulations (symbols) with the single
impurity results. Solid lines are single impurity results which give us the slopes of
the average gallery spacings near :1: = 1. We have used three different values of
KT/ K B as shown in the figure with other parameters fixed; 7A / K A = 73 / K B = 2.0
and KA/KB = 0.1.

 

44

In the EMT for harmonic systems, an effective medium is constructed with
uniform springs of force constant K and equilibrium length he for each intercalant
and K T for the host layers. To study the local distortions caused by an external force
F, the total system is represented by one spring with an effective force constant K c
and an effective equilibrium length he. The effective force constant Ke can be easily

shown to be given by

K
K. _ W’ (3.17)

where the Watson integral W given in Eq. (3.16) contains the combined effects of

all the springs in the system.

To study the effect of randomly intercalated atoms, the theory proceeds by re-
placing one of the effective medium springs, K, by a spring for either A or B
impurity. This system can be visualized as two springs in parallel, one with spring
constant Ke — K and the other with either K A or K 3. The local height and the
energy of either A or B intercalant are obtained by minimizing the energy for each
configuration. With these local heights and the corresponding local distortion en-
ergies, the total energy is expressed as a sum of the energies of A and B impurities
with probabilities 1 — a: and :6, respectively. This total energy is now a function of
two unknown parameters, h,3 and K. Finally, the problem is solved by imposing
two conditions; one is the self-consistent equation for the displacement under the
external force F to obtain K ,, hence K. The other is the minimization of the total

energy to get he.

For the anharmonic system, we similarly consider an effective medium composed
of uniform anharmonic springs characterized by K and 7 with natural lengths h, for

the host-intercalant interaction and harmonic springs characterized by KT for the

45

host layers. To obtain the effective spring parameters Ke and 7,, we consider the

system in the presence of a local external force F. If the external force is applied to

a single site at the origin, the energy can be written as
1 2 1 3
E = EZKM“ — he) — 5270,, - he)

1
+5 ZKTUIn - hm)2 - F(ho - he)-

(3.18)

Minimizing this energy with respect to h,, and performing the Fourier transform on

A,, (E h, — he), we get

=‘I’—(10)1L 2((9) ’

where
<I>(p) E K + 2KT(1— cosp).

The local displacement at the origin A0 can be written as

A0 = 245(1))
= 24%“) +12:A(9)A(P— Q)/¢(P)

= ngfifiZZN q)A(p- q)

FW__(_K) +71A2
K +2K’

1 / K’ is a quantity related to 1/<I>(p) and A(p) through the equation

22,. Z, A(q)A(p - q)/<I>(p)
Z), 23., A(q)A(p - q) '

 

1/K’=

(3.19)

(3.20)

(3.21)

(3.22)

(3.23)

(3.24)

(3.25)

 

46

Different approximation to obtain 1 / K’ will be discussed later. Eq. (3.24) can now
be solved to give the equilibrium length at the origin under the external force F,

i.e.,

 

K’ — \fK” — 27WFK’/K
7 .

 

ho = h. + A0 = h. + (3.26)

Now the same system can also be described by a single anharmonic spring char-
acterized by K,, 7,, and an effective length he under the same external force F.

Then the energy is given by

1 1
E.) = -2—!K.(h — h.)2 — 337.0. — h.)3 — F(h — 1.). (3.27)

Minimizing this energy gives the equilibrium length, h., of the effective single spring

in the presence of F,

Ke — ,/K3 - 2 .F
11, = h, + 7 . (3.28)

‘7:

The parameters for the single effective spring can now be found by equating ha of
Eq. (3.26) with h, of Eq. (3.28) and expanding the resulting equation in terms of

an arbitrary force F. This leads to

K
Kc — W (3.29)
and
7K
7c = W » (3.30)

Like the EMT for harmonic springs, we now replace one of the uniform springs

with K and 7 by either A or B spring. Then the system can be represented by two

47

anharmonic springs in parallel; one with K g = K,3 - K and 7}, = 7e - 7 formed by
removing the spring K and 7 and the other with K a and 7a where a can be either
A or B with probability 1— z and 3:, respectively. The total energy (per site) 8,, for

either of these two types of springs in the presence of F is given by
K’ 7;
£0, = Eflh — he)2 - 5(h — h¢)3 — F(h — he)
K01 7a
+§Ul - hi3.)2 " 370i — h2)3- (3'31)
This energy is minimized with respect to h and gives the local height (ha);

12h. + whifiKé + K...
1!. + 1..

 

(M =

 

_\/(7£he + 73.112 + K; + Ka)’ - ('72 + 7..)(1éhfi + W122 + 2(Kéhe + Kah‘; + F))
7: + 1.. '

 

(3.32)

The average height (h) is therefore

(’1) = (1 — $)<h.4) + Aha). (3.33)

Self-consistency requires that this (h) be the same as h, in Eq. (3.28). The expansion
of the equation h, = (h) up to the second order in F gives three equations for K,

7, and he, namely

 

 

 

l—z :3
he = 7; +7A(9(z‘1)--)/5’TA)+ 7; +7B(.r1(13)-\/7719) _ (3.34)
i _ 1:5,; 335
K. ’ m m <- >
i ._ ___1_$ I x I

 

48

where

9(4) = 72’» + '1th + K2 + K,
f(A) = 731: + 7,119,“ + 2th, + 2K, 119,

TA = 9(4)2 - (72 + 1.4)f(A)
and similar equations with A replaced by B.

Eqs. (3.34)~(3.36) are the EMT equations for the anharmonic model and are
solved numerically in conjunction with Eqs. (3.29) and (3.30) where K, and 7, are
given in terms of K and 7. The results will be discussed in the next section. In the
Appendix B, we show how these three equations recover the earlier EMT equations

obtained in the harmonic system when 7,, = 73 = 0.

Finally we would like to discuss how one obtains l/K’ defined in Eq. (3.25).
Solutions of the above three equations also depend on how we choose 1 / K’ which
appears in the solution to the local distortion problem. As an approximation, we
replace 1/<I>(p) in the second term in the right hand side of Eq. (3.22) by its average
over the p—space which is W/ K , because the most natural way of including the layer
rigidity is to use the Watson integral. Nevertheless, one has to be cautious because
this choice implies that 7, = 7 and therefore 7; = 0 i.e. if we apply a force to
one of the springs of the effective medium (K ,7), the effect of coupling through
KT to the rest of the springs only modifies the harmonic part K —1 K, but not
the anharmonic part. However, it is worth noticing that, regardless of the choice of
1 / K’ , Eqs. (3.34)-(3.36) do recover the results for known cases i.e. for the perfectly
floppy layers (Vegard’s law), for the perfectly rigid layers (Eq. (3.11)) [5], and in the

harmonic potential limit [14] (7 = 0).

 

49

3.5 Numerical Simulations and Comparison
with EMT

Numerical simulations have been carried out for a linear chain with N(= 2000) atoms
for both the LJ and the cubic potential models. Natural length of an atom A(B),
ha(h%), is chosen to be 3.074 (3.820). To generate a configuration of the A1_,B,L
bilayer system, the linear chain is initially filled with all A atoms. Then, A atoms
are randomly selected and replaced by B atoms until the total number of B atoms
N3 = 3N. For a given concentration at, 2000 different configurations are considered.
For each configuration, the conjugate gradient method [20] is used to minimize the

deformation energy given by Eq. (3.1).

Figure 3.4 gives a typical result from the simulations as a: varies from 0 to 1. The
symbols are simulation results and the lines are guide to the eye. We have chosen
K ,4/ K 3 = 0.1, i.e., the smaller atom is more compressible. The middle curve shows

the normalized average gallery spacing of the chain ((1) which is defined by

_ <h> - h°
(d) — W (3.37)

The lower (upper) curve shows the normalized average spacing (dA)((d3)) of the
A(B) atoms. (614) changes much more rapidly with :0 than ((13) because the A

atoms with smaller spring constant cost less energy to compress or expand.

Figure 3.5 compares the numerical simulation results with EMT calculations
for different values of the layer rigidity parameter KT. The solid line is the usual
Vegard’s law given here for the sake of a reference, which is obtained when K T = 0
(perfectly floppy system). The overall agreement between the simulations and the
EMT results is very good. In the region :r —> 1 where small intercalants can be

treated as impurities, EMT results fit extremely well with the simulation results

 

50

 

.IIIIIIIII'IIII'IIII'IIII

 

 

 

 

1.0 '- _a..—(3—EI-
‘ 3-9‘9’8'43 1
C .
0.8;- ._..o “ _
; -"" 1
A 006 _ _
U 1. c-
v ; I
0.4 -— _.
l I
0.2 — -
)- ..
O 0 :1 l l l l l I l l l 1 l l l I l l l 1144 I ll

'00 0.2 0.4 0.6 0.8 1.0
K

Figure 3.4: Average gallery spacing (2:) of a bilayer system, B for B-atom, o for
A-atom in KA/KB = KT/KB = 0.1 and 7A/KA = 7B/KB = 2. The lines are guides
to the eye.

51

 

       
  

 

 

 

:I I I I l I IMTT I I Ifir‘l I I I I I I I I I
1.0 r 1
E x 1
0.6 — f 13 ..
1- , I .,
I I I -(
C l 0’ 1
A 0.6 - 1 ,’ 7
U P I, I 1
V : I I, :
. I , .
0.4 —,'f/ , —
b' I
.; ’ ’ x : x,/x,,=0.01 .
oz '- 0 3 KT/KB=O.1 _..
7 0 = Kr/Ka=1 :
+ : K,/K,=10 *
0.0 1 l l l l l l l l [Ll l l I l l l l l l l l

0.0 0.2 0.4 0.6 0.8 1.0
X

Figure 3.5: Comparison of the numerical simulations (symbols) with the EMT
(dashed curves). The straight solid line represents the Vegard’s law. We have
used four different values of KT / K B as shown in the figure with other parameters

fixed; 7A/KA = 73/KB = 2 and KA/KB = 0.1.

 

I‘.“
'- ‘

 

52

because (1) our choice (7, = 7) handles the anharmonicity of an impurity very
accurately and (2) the anharmonicity of large atoms does not play any significant
role in the gallery expansion as we have discussed in Sec. 3.3. In the small :1: limit,
however, there are many small atoms whose anharmonicity plays important role in
the gallery expansion. Therefore our analytic calculation shows some deviations in
the small a: region due to the several approximations made in obtaining the EMT
results. As seen in the figure, EMT results agree very well with the simulations
for small and large K7 and deviate for intermediate KT values. Deviation between
EMT and numerical simulation results as K T is varied is understood in the following
way. In Eq. (3.25), the quantity 1 / K’ can be thought of as the weighted average
of l/<I>(p) with the weight function w(p) = 2,, E, A(q)A(p - q). We can think of
limiting cases again as in Sec. 3.2. (1) perfectly floppy layers or KT/K = 0,
W(K) = 1: Since <I>(p) is independent of p (see Eq. (3.20)), the weighted average
of 1/<I>(p), 1/K’, is given by l/K regardless of w(p); KT/KB = 0.01 in Figure 3.5
is close to this case. (2) perfectly rigid layers or KT/K —» oo, W(K) = 0: Since
<I>(p) is very large (see Eq. (3.20)), the weighted average of 1/<I>(p), l/K’, can be
considered 0 which is again independent of p. KT/ K 3 = 10 in the figure can be
in this category. Therefore, for these two limiting cases, 1/ K’ is treated accurately

and our EMT results fit numerical simulation results very well for arbitrary values

01 KT.

Figure 3.6 shows gallery expansion as anharmonic parameters are increased for a
given value of K7. All the results of EMT fit well with numerical simulation results
for weak anharmonicity and the agreement is particularly good in the region of a: -—1
1. The EMT results, however, deviate appreciably with increasing anharmonicity

because of the assumptions we have made.

3

53

 

 

 

 

j I I l I I I I I I III I 1 I I I I I I Ii I 1 I I:
1.0 — .. ##gfl-‘I '—

1- In v a
. ,U/er 12:26 j
t ,z" X :
0.8 :- ,:°,)( j
. 1:1,”! X .
P I d
. fl .
A 0.6 — ox’ —
'U ” '1 ‘
V I 17 :
- E”I X : K — =1.0 ‘
0.4 _-_ q 7A/ AWE/KB 1
. 0F 0 : 7A/K‘=73/K3=2.0 .
. 0 2
0.2 L) U 3 7A/KA=73/KB=2'5 __
P, d
'I .
1 I
0.0‘1 l i LLI J l l l l l l J l l l l l l l l I I l l l

0.0 0.2 0.4 0.6 0.8 1.0
X

Figure 3.6: Comparison of the numerical simulations with the EMT results for
different strengths of the anharmonicity. The results from the simulations are de-
noted by different symbols, while those from the EMT are drawn by dashed lines.
KA/KB = 0.1,KT/KB =1.0.

 

 

 

54

3.6 Comparison of Anharmonic Potential Model
with Catchment Area Model

As we have discussed in Sec. 3.1, the normalized gallery height ((1) is given by
(d) = 1 — (1 — (I?)p (3.38)

in the CAM [11, 12], which is an example of an extremely anharmonic (hard sphere)
potential model. In this model, the average gallery spacing evaluated over the small

intercalant sites is
((1,4) = 1 — (1 — z)”l. (3.39)

This is easily obtained from the following argument: the CAM involves hard sphere
atoms, therefore, the large intercalants cannot be compressed. Namely the normal-
ized gallery spacing of the large intercalants, (d3), is always 1 and there is an exact

relation [12] between (d), (01.4), and (d3) given by
(d) = (1 — :r)(d.4) + .7:(d3). (3.40)

From this equation and noting that ((13) = 1, we immediately obtain ((1‘) = 1— (1 —
2:)7". To see how a continuous anharmonic potential describes the z-dependence of
(d) and (d4) as given by the CAM, we fit numerical simulation results obtained with
LJ potential for these quantities to 1 — (1 — .3)" and 1 — (1 — :I:)”'1 (see Figure 3.7).
These four figures are for different values of K A/ K B and for a given value of the
transverse layer rigidity. For a fixed layer rigidity KT/KB, as the ratio K .4 /K 3
decreases, the value of p increases and the agreement between simulations and CAM
results becomes better because expansion of the small atom costs very small energy

similar to the hard sphere case. Note that in the CAM it does not cost any energy

 

55

 

 

 

 

 

 

 

 

1 1
a D f o. o a o o
.. x :6 0 cl P x“ 6 '4
x .o , .
.5 r?
_ X..° II b ’0. q
9'
" l0: '1 l- .“ «4
I? I_.
Hf -- p-2.856 . f -- p-s.429 ‘
t p r l
- ...... p -1 87 l ...... p -2 432
0+ l L l l l l l 1
0 1 0 1

Length

 

 

 

 

 

 

 

 

1 I U . 1
’ XS" 4 L {:8 a
’o' .
r- i; . f- l: u
I? i
.3, '1 .5; .
f —- p-4.078 . -- p-a.548 ,
I
...... '-3. 71 '-5.517
or P. °. of». P.
0 1 0 1

(e) (d)

Composition 1:

Figure 3.7: Comparison of the simulation results obtained with the LJ potential
model to the catchment area model, ((1) = 1 — (1 - 3)”. Four different values of

KA/KB are used [(a) 0.1, (b) 0.07, (c) 0.05, and (d) 0.02]; KT/KB = 0.1.

 

56

to expand the small atoms. Also we find that for a fixed K A / K B, p increases with
the transverse layer rigidity KT / K 3. Thus, the parameter p somehow incorporates
the effects of transverse layer rigidity of the host layers and different intercalant

compressibilities.

In Figure 3.7, we find that the mean gallery separation (d) obtained in the LJ
potential model (simulation results) fits the curve 1 - (1 - 2:)” well even if the
simulation results are slightly below the CAM curve for values of :1: > 0.5. Similarly,
the simulation results for (dA) fit rather well the analytic form 1 — (l — :l:)""'1 for
a: < 0.5 and there are some deviations for a: > 0.5. There is also deviation for
(d3) in the small a: limit. Several conclusions can be drawn from Figure 3.7. First,
each figure shows smaller deviation in the total average gallery spacing than an
average gallery spacing associated with each intercalant (A or B). This can be
explained by the probability of contributions of average spacing of each intercalant
to (d). Namely, as a: —& 1(0), the contribution of small (large) intercalants to the
total average gallery spacing becomes small. Secondly, the gallery expansion of our
system is smaller than that for the hard sphere model because expansion of the
small intercalants to h% costs finite amount of energy unlike the hard sphere model.
Therefore the height of an A-atom is somewhere between 119, and h%, not exactly 11%.
In real systems, it is possible for the intercalants to expand or compress within the
healing length because the intercalants have finite compressibilities. We believe that
if experiments can be performed for the average gallery spacing of each intercalant,
the results will be closer to that obtained via the anharmonic potential model. The
CAM somehow incorporates two important physical effects, the compressibility ratio
of the two atoms and the transverse layer rigidity. However, it is strictly valid in the

limit K A /K 3:0. For finite K A /K 3, (d3) deviates from 1 and ((1,1) is less than 1 in

57

the limit a: —v 1 and these deviations become less significant as the layers become

more rigid.

There is an improved form of the CAM to be denoted as the generalized catch-
ment area model [12]. Since the expansion of atoms costs energy, the equilibrium
length of a small atom lies somewhere between ha and 11% which can be considered

by introducing an additional parameter, a, for ((1,4) besides p. The parameter a is

the value of (d4) in the limit 1: -> 1.
(3,) = 0(1 — (1 — 3)”), . (3.41)

where p’ = p -— 1. The effect of the softening of hard sphere model for a large
atom has been considered. A single large impurity (B atom) has height 19 and is
surrounded by a soft catchment area of small atoms with area p’ = p—l with height
a. If a large atom is surrounded by one or more large atoms, its normalized height

is 1. Therefore,
(d3) = 1 — (1 — mu — z)". (3.42)
From the relation (11) = (1 — 2:)(dA) + w(dg), we get
(d) = 0(1 — (1 — 3)?) + (1 — a) 4 :c — (1 — mm — at)”. (3.43)

Figure 3.8 shows an excellent agreement with numerical simulations of LJ model.
Let us compare the parameters p’ of CAM with that of the generalized CAM. The
generalized CAM allows for A (8) atoms inside the catchment area to have interme-
diate heights 0: (fl). Therefore, within this model, relaxation of A atoms propagates
to a larger distance than that in the CAM. This explains the larger values of p’

(therefore p) in the modified CAM compared to those obtained in the CAM.

The generalized CAM obeys the following relation which can be obtained from

Eqs. (3.41) and (3.42).

58

 

 

 

 

 

 

 

 

 

Length

 

 

    

 

 

1353.421

1 1 l J

(C) (d)

Composition 1:

 

 

 

 

Figure 3.8: Same as Figure 3.7 except that simulation results are fitted with
Eq. (3.41) and Eq. (3.42). The dimensionless catchment area (p’ = p — 1) is in
the figure. Parameters (0,;3) for (d4) and (d3) are as follows. (a) (0.922, 0.937),
(b) (0.940, 0.944), (c) (0.952, 0.950), and (d) (0.973, 0.964).

59

mm) = 113(2) — 113(0)
4.40) 1 - d19(0) .

 

(3.44)

Figure 3.9 is a plot of this relation with LJ simulation data. As K A / K 3 decreases,
the LJ results show deviation from the relation because the generalized CAM treats
the softness of large and small atoms in a symmetric way whereas in the case of LJ
simulation, the asymmetric shape of anharmonic potential gives different amount of
expansion and contraction. This suggests that the effect of the softness of the large

and small atoms should be considered differently. As K A / K 3 decreases, this effect

is more pronounced as we see from the figure.

Although we have compared the generalized CAM with LJ simulation keeping
KT/KB fixed and varying K A / K 3, we believe similar agreement between the two
will be obtained for different KT/ K 3 but fixed K A / K 3.

3.7 Conclusion

In summary, we have studied the gallery expansion of ternary intercalated systems
within anharmonic potential models. We find that the anharmonicity is important
in layered ternary systems, especially for low concentration of small intercalants (i.e.
:r -+ 1). When the compressibilities of the two intercalants are the same, one obtains
Vegard’s law within a harmonic model. On the other hand, the inclusion of anhar-
monicity gives rise to deviations from Vegard’s law. These deviations increase with
increase in the transverse layer rigidity. Effective medium approximation appears
to be adequate for small anharmonicity, and in the limit of small and large layer
rigidity. The results of the “Catchment Area Model” can be understood within a
Lennard Jones model and simple generalization to the CAM can be made to fit the

gallery expansion obtained within LJ model very well.

 

60

 

 

 

 

 

 

 

 

 

 

 

 

.- 88 a .- 0 ..
. 3 . .. 9 .
. 3 . . 9 .
.. - P. -1
of- l J 1 1.1 or I l l l
A 0 1 0 l
4;” (a) (b)
E
lr—‘l’ I 1 a 1—1' I 18;“
898 fl 0388
.- O .. .. on ..
h a d I- a ..
O
I- II To II
a O
0T— l L 1 l o?— l l 1 J

 

 

 

 

(e) (d)

Composition it

Figure 3.9: Plot of Eq. (3.44) by using LJ simulation results of (dA) and (d3).

 

Bibliography

[1] F. S. Galasso, Structure and Properties of Inorganic Solids (Pergammon, New
York, 1970).

[2] L. Vegard, Z. Phys. 5, 17 (1921).
[3] Y. Cai and M.F. Thorpe, Phys. Rev. B 46, 15879 (1992).

[4] J. B. Thompson, in Researches in Geochemistry, Vol. II, ed. by P. H. Abelson,
(wiley, New York, 1967) p. 340.

[5] J. R. Dahn, D. C. Dahn, and R. R. Haering, Solid State Comm. 42, 179 (1982).
[6] J. E. Fisher and H. J. Kim, Phys. Rev. B 35, 3295 (1987).

[7] S. A. Solin and H. Zabel, Adv..Phys. 37, 87 (1988).

[8] P. C. Chow and H. Zabel, Phys. Rev. B 38, 12837 (1988).

[9] D. Medjahed, R. Merlin, and R. Clarks, Phys. Rev. B 36, 9345 (1987).

[10] W. Jin and S. D. Mahanti, Phys. Rev. B 37, 8647 (1988).

[11] H. Kim, W. Jin, S. Lee, P. Zhou, T. J. Pinnavaia, S. D. Mahanti, and S. A.
Solin, Phys. Rev. Lett. 60, 2168 (1988); S. Lee, S. A. Solin, W. Jin and S. D.

Mahanti, in Graphite Intercalation Compounds: Science and Technology, Ed. by

61

 

11

62

M. Endoh, M. S. Dressellhaus and G. Dressellhaus, p. 41 (Materials Research
Society, 1988).

[12] M. F. Thorpe, Phys. Rev. B 39, 10370 (1989); unpublished.

[13] S. Lee, H. Miyazaki, S. D. Mahanti, and S. A. Solin, Phy. Rev. Lett. 62, 3066
(1989).

[14] M.F. Thorpe, W. Jin, and SD. Mahanti, Phys. Rev. B 40, 10294 (1989); Dis-
order in Condensed Matter Physics, Ed. by J .A. Blackman and J. Tagiiefia, p.22
(Clarendon, Oxford, 1991); M.F. Thorpe and SD. Mahanti, Lectures given at
NATO Summer School, “Chemical Physics of Intercalation II” , Chateau de Bonas,
France, 1992; Lectures given at the 7th International School on Condensed Matter

Physics, Varna, Bulgaria, 1992.

[15] Y. Cai, J. S. Chung, M. F. Thorpe and S. D. Mahanti, Phys. Rev. B 42, 8827
(1990).

[16] Hyangsuk Seong, Jean S. Chung, and S. D. Mahanti, Bull. Am. Phys. Soc. 35,
781 (1990); preprint.

[17] N. Bernardes, Phys. Rev. 112, 1534 (1958).

[18] C. Kittel, Introduction to solid state physics, (5th ed., McGraw-Hill, 1977)
p. 77 .

[19] S. Feng, M. F. Thorpe, and E. Garboczi, Phys. Rev. B 31, 276 (1985).

[20] W. M. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vettering, Numerical
Recipes (Cambridge University Press, 1986) p. 301.

[21] T. Horiguchi, J. Math. Phys. 13, 1411 (1972).

 

Chapter 4

Structure of 2-Dimensional (2D)
Repulsive Screened Coulomb
Systems / Model Intercalation
Compounds

4.1 Introduction

Structural and dynamic properties of systems with competing interactions continue
to be one of the most interesting areas of condensed matter physics [1]. Of particular
interest over the past decade has been the study of structural properties of binary
and ternary graphite intercalation compounds (GICs) [2, 3]. The interplay be
tween intercalant-intercalant interactions and graphite-intercalant interactions are
ultimately responsible for the intercalation process. Furthermore, depending on
these competing interactions which are naturally sensitive to the choice of inter-
calants and their planar densities, various structural features appear. Considerable.
amount of experimental [2, 3] and theoretical [4] works have helped in elucidat-
ing many of the interesting physical properties of these systems. The stage n 2 2
GICs have provided rich sources for the study of structural properties and phase

transitions in two dimensions. However, one area where a detailed microscopic

63

64

understanding is still lacking is the incommensurability of stage-n(n Z 2) GICs
containing Cs, Rb, and K atoms as intercalants. We address this issue by probing
the structure of domains and domain walls using Molecular Dynamics (MD) simu-
lations and then propose a simple model for the low temperature structure of the
alkali GICs. Although both X-ray and neutron diffraction studies have been cru-
cial to the elucidation of the structure, in this chapter we will consider only X-ray
diffraction experiments since both yield the same results on structures. In addition,
quasi-elastic and inelastic neutron scattering studies have been of great importance

for understanding the dynamics of (3103 and will be discussed in chapters 6 and 7.

4.2 Physical Systems

Pure graphite is characterized by a layered structure, originating from a strong co-
valent bonding of the carbon atoms in the plane and weak Van der Waals coupling
between the planes. The in-plane C-C bond length is 1.42A, and the interlayer
spacing is 3.35/4 (see the upper panel in Figure 4.1). The planes are arranged to-
gether in the sequence ABA and are related to each other by the translation vector
g = (2a - b)/3, where a and b are basis vectors of the graphite unit cell. Un-
der appropriate conditions of graphite temperature and alkali vapor pressure [5],
intercalation of alkali atoms into the graphite interlayer spacings occurs. The in-
tercalated atoms lead to three main structural changes of the graphite hosts. First,
‘staging’ where the alkali atoms intercalate in a regular fashion so that an inter-
calated layer is followed by a constant number (n) of graphite planes, i.e. stage-n
GICs are characterized by a stacking sequence of graphitic and intercalant layers in
which the neighboring intercalant layers are separated by n graphitic layers (see the

lower panel in Figure 4.1). Second, alkali layers are always flanked by equivalent car-

65

 

 

 

 

'9' arm m sun in arm on sun

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

O... C... C...‘ O...‘

‘0... . .
0...

case ‘_ ‘

A O

 

0 “HIGH “(3160' LIVE!

0 O O 0 LAYER 0' INYCICCLAYI
“LIAM NIT“. ATOII

Figure 4.1: Upper panel is a structure of graphite. Lower panel shows stacking
sequence of carbon and alkali layers.

66

bon planes, thereby causing a rearrangement of the stacking order of the graphite.
Third, the interlayer spacing of adjacent graphite planes increases considerably (e.g.
from 3.35 A to 5.71 A for Rb) upon intercalation, while the in-plane structure of
the graphite remains almost unaffected. During intercalation, approximately one
electron per alkali atom is transferred to the graphite. Experiments indicate that
these extra electrons not only compensate for the charge of the remaining ionic al-
kali layers but also screen them [6]. The intercalant ions interact predominantly via
screened Coulomb repulsion and feel an effective one particle corrugation potential
produced by the graphite layers. Thus, these systems exhibit a competition between
two length scales, periodicity of the substrate corrugation potential and the aver-
age separation between the intercalant ions controlled by their planar density (p).
In addition, there is also a competition between two energy scales, the interaction
energy/ particle and the depth of the corrugation potential which is determined by

p and the strength of the interparticle interaction.

In stage-1 GICs of K, Rb, and Cs, the intercalants form a commensurate trian-
gular (2 x 2)R0° structure where R0° and (2 x 2) imply that the unit cell vectors
of the intercalant superlattice are parallel to the graphite unit cell (a, b) and are
twice as large. The layer stoichiometry is A03, where A stands for the intercalant
and C for the carbon (note there are two carbon atoms per graphite unit cell). This
gives p = 1/4 (in units of 1 / area of the graphite unit cell). In contrast, structures of
stage-n (n 2 2) systems are more intriguing. Chemical stoichiometry of these com-
pounds is in general A012,, and the layer stoichiometry is A012, giving p = 1/6. In
the absence of the corrugation potential, the ground state is a triangular lattice with
an arbitrary orientational epitaxy angle ¢ (measured with respect to graphite recip-

rocal lattice vector a'-axis, also known as Novaco-McTague angle [7]). A question

67

of fundamental interest is the effect of corrugation both on the intercalant structure

and the intercalant dynamics.

The liquid state structure in stage-2 compounds is well understood through care-
ful MD work of Moss et al. [8, 9, 10] and Chen et al. [11] in terms of a highly
structured liquid close to the liquid-solid phase transition. The nature of the solid
structure is, however, not as well understood although qualitative models have been
proposed by Clarke et al. [12, 13] for A024(A = Cs, Rb, and K) in terms of different
types of commensurate domains separated by discommensurations (domain walls).
Also Ginzburg-Landau models have been proposed to understand the domain struc-
ture [14]-[17]. From X-ray diffraction measurements in CsC24 and RbCu, Clarke
et al. [12, 13] suggested the coexistence of commensurate (J7 x J7) and (2 x 2)
regions. They explained their diffraction results for Cs by postulating the coexis-
tence of (J7 x J7) R¢ domains with (6 = (+/-—)19.11°. The commensurate domains
were of size ~ 40A containing about 30 Cs atoms along with narrow domain walls
oriented parallel to the graphite (110) directions (see Figure 3 in ref. [12]). They,
however, neither discussed the precise nature of the domain walls nor analyzed the
spatial distributions and relative amounts of (J7 x J7) and (2 x 2) domains, the
latter being determined by the density constraint. In the case of RbC'24, it was not

quite clear whether (2 x 2) domains were present or not.

Following the suggestions of Clarke et al. [12], several theoretical studies based
on Ginzburg-Landau theory of discommensuration in the incommensurate phase
were carried out pointing out the importance of third order terms in stabilizing
the (W x J7) domains [14]-[l7]. These theories were able to explain the observed
Novaco-McTague angle d) and intensities of the higher order superlattice reflections

semiquantitatively, although they did not address the structure of the domain walls

68

in detail. As we will show later, domain and domain wall sizes in these systems
are comparable and one should therefore also know the microscopic structure of

the walls.

4.3 Potential Model

MD simulations were carried out for RbCu using the potential model of Ref. [11].
The one particle corrugation potential in this model was obtained by Moss et al. [8]

by comparing the liquid structure factor S (k) with X-ray experiment and is given

by
V1(r) = —2K[2 cos(21ra:/a) cos(21ry/J3a) + cos(41ry/\/3a)]. (4.1)

with 2K = 0.9 (in units of 30013;) i.e. 2K/k3 = 270K, and a = 2.46.74 which is the
graphite unit cell constant. The two-body potential is given by
V2 = gazeXM-I‘ml/ r15, (4-2)
:1
with q = 4.8028 x 10'10 can and I‘ = 0.4921”.

The above potential has successfully reproduced the liquid structure data at
250K. Our MD simulations using constant energy algorithms agree with the earlier
results of Chen et al. [11]. Before discussing the solid structure results, let us
discuss some energetics. The strength of the corrugation potential as measured by
V,‘“‘“" — VIM" = 1080K i.e. we are in the strong to intermediate corrugation limit for
T < 300K. The energy/ particle for different commensurate structures are 17742K
(J3 x J3); 9825K (2 x 2); 4476K (2 x 3); 3114K (J? x J7). Although the (2 x 3)
structure has the correct layer density, the lower energy of the (J7 x J7) structure

favours the latter.

69

We can generalize this problem by changing parameters which appear at
Eqs. (4.1) and (4.2). When 2K = 0, the system corresponds to a purely 2D re-
pulsive screened Coulomb system. As 2K is increased, the effect of host layers (or
substrates) increases and how the physical properties change with 2K will be dis-
cussed in chapters 5-7. As we change I‘, the interaction between intercalants (or
adsorbates) varies from purely Coulomb (I‘ = 0, long range) interaction to highly
screened Coulomb (I‘ -r co, very short range) interaction. This chapter will focus

on a system with 2K = 0.9 and I‘ = 0.49 appropriate for stage-2 Rb GIC [11].

4.4 Details of Molecular Dynamics Simulations

MD simulations were performed with 216 and with 864 Rb ions distributed to insure
a planar density of p = 0.0318 Rb ions/A2 for the stage-2 Rng. system. In addition
to this planar density, we have also studied a 301 particle system with a slightly
smaller density of p (0.0311 Rb ions/A2) which corresponds to the stage—2 RbC24,57
system to see the effect of change in the density on the ground state structure. MD
simulations were carried out using a sixth order Gear predictorcorrector algorithm
[16] to solve the Newton’s equations with periodic boundary conditions. The time
step was taken to be 0.0029 picosecond (ps) and typical equilibration times of at
least 350 ps were used for 216- and 301-particle systems and 1 nanosecond (us) for
864-particle systems, respectively. The energy/ion was noted to be conserved to
99.994% accuracy over 150 ps run for these systems. The range of the two-body
potential was kept within a circle of radius 27.06A centered on any Rb ions. The

results remained essentially the same when a larger cut off radius was used.

Most of the simulations were carried out in the temperature range 100K < T <

300K and in some cases we cooled down the system to as low as 3K. The 216-particle

70

system was initially heated from a triangular structure and allowed to melt and to
equilibrate. It was then slowly cooled in temperature intervals of 10K through
the transition region. Finally the cooled solid was slowly heated to see whether
liquid-solid transition showed hysteresis or not. For the RbCz4,5-, system we chose
a different initial configuration in order not to bias our simulations towards any
particular ordered structure as its ground state. For this system we started from a
random configuration; this corresponded to an extremely hot nonequilibrium system.
The system was cooled down rather quickly to 300K from the initial nonequilibrium
temperature of ~ 10‘K and was allowed to equilibrate at that temperature for
sufficiently long time. We then slowly cooled the system by temperature intervals
of 10K until the ground state was reached. It turned out to be a Periodic Domain
Wall solid. This solid was then slowly heated to probe the solid-liquid transition.
For these systems, all results of cooling and heating runs were compared to ensure

that thermal equilibrium was achieved.

Here we report our results for two different temperatures, 250K (liquid) and 3K
(solid) [18]. The study of the melting transition will be given in chapter 5 [19] and

dynamics of these systems will be studied in chapters 6 and 7.

4.5 Structural Properties

The structure of simple monatomic fluids is characterized by a set of distribution
functions for the atomic positions, the simplest of which is the pair distribution
function g(r) which gives the probability of finding a pair of atoms a distance r

apart.

g(r) = 17.42333 — n.» (4.3)

 

71

where A is the area of the MD cell, N is the number of Rb ions in the cell and
the angular brackets imply a time average. In Figure 4.2 (a) and (b) we give the
spherically averaged real space pair distribution function g(r) for the solid and the
liquid, respectively. Clearly, the liquid is strongly perturbed by the host layer at
this temperature and is structured, as noted earlier by Fan et al. [9] and Chen et

al. [11] following the earlier suggestion by Parry [20] for 03024.

The solid phase g(r) shows three distinct peaks at 2.19, 2.53, and 2.97 in units V
of graphite unit cell length and these can be related to distances of 2, J7, and 3,
respectively [21]. A similar three-peak structure was obtained from a Monte Carlo
study by Plischke et al. [22] and suggests that the system may show coexisting
domains of (J7 x J7) and (2 x 2) consistent with the picture proposed by Clarke et
al. [12] for 03024. The positions of the intercalants joined by lines connecting near-
est neighbours using Voronoi polyhedra constructions [23] are shown in Figure 4.3
(a) and (b) corresponding to two different concentrations of Rb atoms (Rng4 and
Rng4,57). In contrast to the picture of Clarke et al. for 03024 (presence of large
domains of (J7 x J7) structure containing about 30 atoms oriented in different
directions and separated by narrow domain walls of unknown atomic structure),
we see placquettes of (J7 x J7) structure consisting mostly of 7 atoms (we will
call these nano-domains) for Rng4 (Figure 4.3 (a)). The intersection of three such
placquettes forms a triangular placquette of (2’ x 2) structure containing 3 atoms.
The walls between two (J7 x J7) nano—domains consist of (2 x 3 x J7) triangular
placquettes which give rise to the peak near 3 in g(r). These domain walls should
also give a peak corresponding to r near J1—3 = 3.6 which we also find in g(r) (see
Figure 4.2(b)). Our MD simulations give a picture rather similar to the one sug-

gested by Zabel et al. [24] for Rng4 (see their Figure 1 and figures in ref. [25]).

 

72

 

 

 

 

 

 

 

 

 

 

12
" 0
I (a)
10—
.L [ .
C
g(r) 6 L— [T
t .. "
4- ".. l. .. l ["1

 

 

 

 

 

Figure 4.2: Pair correlation function g(r) for RbC24 as a function of r (in units of
graphite lattice constant a=2.46A) (a) for the solid at 3K (upper panel) for the 216
particle system and (b) for the liquid at 250K (lower panel). Note that the zero for
the solid is 3.

73

 

 

 

)

(b

 

 

a) R6024 and (b) Rng4_57. The

(

solid, dotted, and dashed lines correspond to the nearest neighbor distances near

Figure 4.3: Low temperature solid structure of

7,2, and 3 (in units of graphite lattice constant), respectively.

74

In fact, by choosing a slightly smaller density of Rb ions we see a perfectly ordered
arrangement of J7 x J7 domains with (2 x ’3) domain walls.

To compare the MD results with X-ray diffraction measurements we calculate
the structure factor S (k) for the atomic positions given in Figure 4.3 (a). The static

structure factor is defined as follows,

2) , (4.4)

where rj(t) is the position of the j-th particle at time t. The structure factor for

2 ask-1’50)

J

 

 

S(k)= jlv-<

the solid phase (216 particle system) is given in Figure 4.4. Our MD simulation
shows two prominent peaks (I; and [2) at the same magnitude of wave vector k, =
1.21A’1 but different angles with respect to the graphite reciprocal lattice vector
a“-axis. The angle associated with the dominant peak 11, denoted as ()3, is 102°
whereas that associated with the weak peak I; which we denote as :12 is 35.8°. The
strengths of these two peaks are 0.52 and 0.052, respectively. There is a third peak
13 (k = 1.4021“, angle with the a‘-axis is 7.6°) which is very weak (strength 0.015)
and will not be discussed any further. The 1:, values and angles of the two strong
peaks measured from the graphite a‘-axis are given in Table 4.1 along with the
available experimental values and predictions of the Periodic Domain Wall Model

to be described in the next section.

Clarke et al. [12] saw two peaks similar to 11 and 12 in 03024 (1 and 2 in
their notation in ref. [12]). The corresponding 1:, and angles 43 and 1b are given
in the Table 4.1. In addition, they also saw a peak as strong as I; (1’ in their
notation) which originated from the rotational equivalent of 11 corresponding to the
Novaco McTague angle —¢. The origin of the low intensity peak I; was ascribed by

Clarke et al. to the inter-domain scattering associated with the coexistence of two

75

 

 

 

 

 

 

 

Figure 4.4: Structure factor S (k) for k in the let-k, plane for RbCu at 3K obtained
from MD simulations for the 216 particle system. Only the dominant peaks 1;, I2
(see Ref. [18]) and their six-fold symmetry counterparts are clearly visible.

76

Table 4.1: Comparison of angles and peak positions of S (k) between Experiment
(EXP), Molecular Dynamics (MD), and Periodic Domain Wall Model (PDW) (see
text). 45 = N ovaco McTague angle; ¢ = angle of the peak produced by inter-domain
scattering (angle of I; in Ref. [13]); 1:, is the magnitude of k at angles 4: and (1:.

¢ (¢) k»
EXP MD PDW EXP MD PDW

R5024 10.1 (35.5) 10.2 (35.8) 11.5(33.3)‘ 1.22 1.21 1.19.
6302, 14.5 (28) ? (?) 14.5(27.6)b 1.16 ? 1.15“

 

 

 

 

 

‘ for L=M=1 (RbC-um)

b for L=2, M=1 (Cngm)

different orientations of (J7 x J7) commensurate domains. Our MD simulations
show only one orientation of these domains (either 1 or 1') thus indicating that it
is not necessary to invoke domains with different orientations to explain the origin
of the I: peak. It is sufficient to have domain walls, in our case the two domains
sandwiching a wall being simply shifted with respect to one another by two units of
a graphite lattice vector. For RbCu the dominant peak 11(lc, = 1.22A",¢ = 101°)
[12] is in good agreement with our MD results. The lower intensity peak (12) for
this system was not discussed in Ref. [13] but the corresponding 1:, and the angle 1b

agree very well with the X-ray diffraction results of Rousseaux et al. [27].

The liquid state S (k) at 250K (in Figure 4.5) shows clearly the effect of graphite
host layer with the dominant anisotropic peak at k = 1.2:?1 corresponding to a
corrugated liquid [26]. The inset in Figure 4.5 shows the S (k) averaged over angle.

The sharp peak near 3A"1 is a result of the graphite corrugation potential.

77

 

 

 

_....l....,r...,:
32— —I

S(k)2:— —
13— f
05 ..1.......I.EI
0 2

4
k(A"1)

   

4 o '0
.0113...»
' I in '

   

  

. a
- 4"
o -- l A
o u o l
0.,”20’!’
0...:

 

 

 

Figure 4.5: The liquid structure factor 5' (k) for k in the kx-ky plane for RbCu near
T=250K shows an anisotropic ridge due to the effect of corrugation. The inset is
S(k) averaged over angle.

78

4.6 Periodic Domain Wall Model

To understand more about the physical origin of the two peaks discussed in the
previous section, we propose a Periodic Domain Wall (PDW) model which consists
of periodic arrays of commensurate (J7 x J7) domains of total width 2L (all of
same orientation) interspersed by domain walls of width M (in units of 2a). The
domain walls consist of (2 x 3) regions and the region where three domains meet
consists of a triangular array of (2 x 2) structure (see Figure 4.3(b)). For arbitrary
(L, M) the super-lattice unit cell is obtained by joining the centers of the (J7 x J7)

domains and the unit cell vectors are given by

A = 2(2L + M)a + Lb, (4.5)

= —La + (5L + 2M)b, (4.6)

where a = a5: and b = (a/2):i: + (J3a/2)y, :i: being a unit vector along the (10)

direction of a graphite in real space and 9 being a unit vector perpendicular to :i: in

 

a right handed coordinated frame. The size of super-cell is aJ2le + 18ML + 4M3.
One often considers this as the domain size due to unknown nature of the narrow
domain-wall structure. However, in our case, domain is clearly defined and its size is
obtained by putting M=0 i.e. aLJ2T. The reciprocal lattice vectors corresponding

to the super-lattice unit. cell vectors are

 

 

 

_ 21r(5L + 2M) , 3L + 2M ..
K" ’ 31(21L2 + 18LM + 4M2)(‘” " fig}, + 2M)”)’ (4'7)
21rL , 9L + 4M ,
KB —a(21L2 + 18LM + 4M2)” " (5L 3’)“ (4‘8)

The density in the PDW model is determined by (L, M). The number (N31,) of Rb
ions in the super-cell is 3L(L+M)+M2 where as the number (N0) of Carbon atoms

in it is 2(21L2+18ML+4M2). Therefore, we can obtain a planar density from N31,

 

79

and No for given (L, M). Also m = Nc/Nm can be used for the nomenclature of
stage-n 2 2 GICs as AC“. where A stands for the intercalants. In this model one
goes from the commensurate (J7 x J7) structure (M=0, arbitrary L) whose layer
stoichiometry is AC“ to the commensurate (2 x 2) structure (arbitrary M, L=0)

corresponding to AC3.

Our MD simulations suggest that RbCu can be described very well by adding a
small number of Rb atoms to the PDW model with L=M=1 (this corresponds to the
stoichiometry RbCusr). To check this idea, we studied a 301 particle system with
density of Rb ions corresponding to Rng4,57. Figure 4.3(b) gives its low temperature
structure. The structure factor for the defect-free case has been calculated exactly
and the corresponding values of hp, angles, and the strengths of the two peaks [1 and
I; are respectively (1.19, 1.19)A‘1; (11.5°, 33.3"); (0.573, 0.137) (see Table 4.1). The
larger value of the intensity of the I; peak in the PDW model (0.137) for RbC24,57
compared to the MD results (0.052) for RbCu is due to the perfect ordering of the
Rb ions in the PDW model.

We find by actual calculation that d (the Novaco-McTague angle) and 1,!) (the
angle of the peak produced by inter-domain scattering, angle of I: in Ref. [13])

satisfy the relation
sin .1) = fisinwo — is), (4.9)

for the PDW model with M=1 and arbitrary L. In the limit L -» 00, this relation
gives (11 = 1/2 = 19.1° i.e., the two peaks merge with each other for the commensurate
(J7 x J7) case. In Figure 4.6, we plot 1]) and 43 for different values of k, which is
a measure of the layer density and is directly related to L and M by the reciprocal

lattice vectors. This finding confirms that the appearance of the peak 12 originally

 

80

 

 

 

 

40 I I I I I I 17 T I 1M I U I I I U 1 I 1 T I U
I- I D ..
’0? f p '1
0.) 1°
0 ’ I/ r
:0 so - ‘15 /‘ 4
o . .
3 . .fd .
g3) . f -
II-t .. ..
2? ._
a1 20 #
>. F J
x r-
:33 . o o i
a. - 0 D .
m 10 r U i D —
.1 l l l l l l l l l l l l l I l l J 1 ii I 1 ll
0.37 0.38 0.39 0.40 0.41 0.42

kp/a'

Figure 4.6: Novaco McTague angle 43 and the epitaxy angle 1b corresponding to the
peak produced by inter-domain scattering as functions of domain size L with M=1
obtained within periodic domain wall model. Here, It, along the x-axis varies as
L and M change. Solid squares are our MD simulation results for Rng4. White
circles (squares) are various experimental results [2, 15] for Rng4(Cng4). A starlike
symbol is an experimental result for C 3036.

 

81

noted by Clarke et al. [12] is unrelated to the presence of domains of c0mmensurate
(J7 x J7) region with two different orientations, namely :tl9.11°. Indeed, such
configurations which generally require the Rb ions to be located closer than 2a in
the regions between the domain walls (if we stick to the layer density of RbCu)
are energetically unfavorable in view of the strongly repulsive two body potential
between these ions. Furthermore, the 03034 results of Clarke et al. [12] can be
approximated by the L=2, M=1 PDW model which corresponds to a stoichiometry
03026.1, and gives 1: = 1.15.21", ()3 = 14.5°,1/2 = 27.6° in excellent agreements with
the experimental values 1.16 A“,14.5°, and 28° respectively. The corresponding
intensities in the PDW model are respectively 0.568 and 0.144 and the number of

Cs atoms in a domain is 19.

4.7 Summary and Conclusion

In summary, (i) the corrugation potential obtained by comparing the structure factor
with X-ray scattering experiment in the liquid state can be successfully applied
to the study of the low temperature microstructure in alkali GICs. The success
of the chosen potential is confirmed by the agreement between the experimental
S (k) and that obtained from our simulations for Rng4. (ii) We also find that the
simulated S (k) for the stage-2 GIC closely agrees with the experimental S (k) for
stage-3 systems [13]. This observation confirms our original assumption concerning
weak interlayer coupling and effective two-dimensionality of these alkali GICs. (iii)
Simulations suggest that the real space structure of the Rb intercalants consists of
thick domain walls lining the domains with length scales comparable to the domain
walls. (iv) Finally, our S(k) calculations for the L=2, M=1 PDW model which

recognizes the importance of the domain walls shows remarkable agreement with

82

experimental S (k) for the stage—2 Cs GIC. This last success indicates that the PDW
model describes very well the main features of the low temperature structure of the

stage-2 and higher GICs of Rb and Cs.

 

Bibliography

[1] Competing Interactions and Microstructures: Statics and Dynamics, R. LeSar,
A. Bishop, R. Hefner, eds., Springer (1988).

 

[2] Graphite Intercalation Compounds, Vol 1: Structure and Dynamics, H. Zabel

and S. A. Solin, eds. Springer Series on Topics in Current Physics, Springer (1990).
[3] S. A. Solin and H. Zabel, Advances in Physics, 37, 87 (1988).

[4] S. A. Safran, Solid State Physics vol. 40, ed. by H. Ehrenreich, and D. 'Ilurnbull,
Academic Press (NY) (1987).

[5] RC. Eklund, Intercalation in Layered Materials, ed. M.S. Dresselhaus, Plenum,
New York, p163 (1986). l

[6] J.E. Fisher, Physica, 9913, 383, (1980).
[7] A. D. Novaco and J. P. McTague, Phys. Rev. Lett. 38, 1286 (1977).
[8] S. C. Moss, et al., Phys. Rev. Lett. 57, 3191 (1986).

[9] J. D. Fan, Omar A. Karim, G. Reiter, and S. C. Moss, Phys. Rev. B 39, 6111
(1989).

[10] J. D. Fan, George Reiter, and S. C. Moss, Phys. Rev. Lett. 64,188 (1990).

83

84

[ll] Zhuo—Min Chen, Omar A. Karim, and B. Montgomery Pettitt, J .‘Chem. Phys.
89(2), 1042 (1988).

[12] Roy Clarke, J. N. Gray, H. Homma, and M. J. Winokur, Phys. Rev. Lett. 47,
1407 (1981).

[13] M. J. Winokur and R. Clarke, Phys. Rev. Lett. 54, 811 (1985).
[14] Y. Yamada and I. Naiki, Jour. Phys. Soc. of Japan, 51, 2174 (1982).
[15] M. Suzuki and H. Suematsu, Jour. Phys. Soc. of Japan 52, 2761 (1983).

[16] M. P. Allen and D. J. Tildesley, Computer Simulations in Liquids, Clarendone,
Oxford (1987).

[17] M. Suzuki, Phys. Rev. B 33, 1386 (1986).

[18] Hyangsuk Seong, Surajit Sen, Tahir Cagin, and S. D. Mahanti Phys. Rev. B
45, 8841 (1992).

[19] Hyangsuk Seong, S. D. Mahanti, Surajit Sen, and Tahir Cagin, Phys. Rev. B
48 (1992).

[20] G. Parry, Mat. Sci. Eng. 31, 99 (1977).

[21] Surajit Sen, Tahir Cagin, Hyangsuk Seong, and S. D. Mahanti, Recent Devel-
opments in Computer Simulation Studies in Condensed Matter Physics IV, ed.
by D. P. Landau et al., Springer (1991).

[22] M. Plischke and W. D. Leckie, Can. J. Phys. 60, 1139 (1982).
[23] GP. Voronoi, J. Reine, Angew. Math. 134, 198 (1908).

[24] H. Zabel, et al., Phys. Rev. Lett. 57, 2041 (1986).

 

85

[25] S. C. Moss and R. Moret in Ref.2.

[26] H. Zabel, A. Magerl, J. J. Rush, and M. E. Misenheimer, Phys. Rev. B 40,
7616 (1989).

[27] F. Rousseaux, R. Moret, D. Guerard, and P. Lagrange, Phys. B 42, 725 (1990);
F. Rousseaux, et al., Synth. Met. 12, 45 (1985).

 

Chapter 5

Melting of a Repulsive Screened
Coulomb System in 2D - Effect of
Corrugation

5.1' Introduction

Physical properties of two—dimensional (2D) solids and liquids and the nature of the
solid-liquid transition continue to be of theoretical and experimental interest since
the dislocation- and disclination-mediated melting ideas of Kosterlitz and Thouless
[1], Nelson and Halperin [2], and Young [3]-[6]. The nature of the melting transition
(whether first order driven by grain boundary melting [7] or continuous driven by
unbinding of above topological defects) and the existence of the intermediate ‘hexatic
phase’ with quasi long-range positional order and long-range bond orientational
order sandwiched between the ‘2D solid’ phase and the ‘2D isotropic liquid’ phase

are still uncertain [4].

In real physical systems such as atoms adsorbed on 2D substrates or interca-
lation compounds with sufficiently weak interlayer interaction which renders these
systems quasi two-dimensional, the dominant effect of the substrate is the periodic

corrugation potential. The latter affects the structures of both the solid and the

86

 

87

liquid phases and the solid-liquid phase transition [2, 8, 9]. In particular, one has
to take into account both the commensurability of the periodic potential and the
strength of the corrugation. In the limit of weak corrugation, Nelson and Halperin
[2] have argued that for an incommensurate potential the spatial variation of the
Novaco-McTague angle [10] gives rise to a new elastic constant (7) in addition to
the two usual Lamé constants (A and )1) associated with a 2D isotropic elastic solid.
The elastic constant 7 determines the energy related to the spatial variation of the
bond angle field in the presence of a finite substrate potential. Thus the corrugation
potential acts as a field which couples to the bond orientational order parameter.
Consequently, one has only the dislocation unbinding transition and the disclination
unbinding transition is washed out just as a ferromagnetic transition gets washed
out in the presence of a uniform magnetic field. For intermediate to strong incom-
mensurate corrugation potential the precise nature of the transition is not known

(see below the discussion in the lattice-gas limit).

When the corrugation potential is commensurate, the nature of the melt-
ing transition is profoundly affected. For coarse substrate mesh, the Kosterlitz-
Thouless picture becomes inappropriate and even for a weak substrate potential a
lattice-gas picture describes the solid-liquid phase transition [2]. For sufficiently fine
substrate corrugation, the commensurate solid undergoes a transition to a float-
ing solid phase which then melts through dislocation unbinding. In this case the
disclination unbinding transition is also washed out. In the limit of strong corru-
gation potential the system of course behaves like a lattice gas [11]. Even in this
case Nelson and Halperin speculated that one could have a melting process similar
to that for a continuum model if the corrugation mesh size was sufficiently small.

Thus one expects to see very different melting behaviour depending on the strength

 

88

and the periodicity of the corrugation potential. If on the other hand the melting
transition is determined by the grain boundary condensation [7] one would like to
know the effects of the corrugation on the grain boundary formation and hence on

the transition.

Physical systems such as physisorbed rare gas atoms on graphite [6, 12, 13],
H2(D2) on graphite [14], and graphite intercalation compounds [15] fall in the weak
to intermediate corrugation strength limit. In these systems the wavelength of the
corrugation potential, a, is typically r(J3 S r 5 J7) times smaller compared to
the systems periodicity (A0) in the absence of corrugation. We do not know of
any real physical system where a > A0. Thus to understand the full scope of
melting on a corrugated surface for different corrugation strengths and different
substrate periodicities it is necessary to undertake a systematic simulation study of
the melting phenomena for systems with different densities and different strengths

of the corrugation.

Recently, Vives and Lindgard [16] have carried out extensive Monte Carlo simula-
tions on a commensurate (J3 x J3) Lennard-Jones (LJ) system to study the nature
of the transition and the decay of translational and bond orientational correlations
as a function of the strength of corrugation. They have also extended these calcula-
tions to the incommensurate case appropriate for the H2(Dg)/graphite system [17]
although the major emphasis in this work was to understand the low temperature
structure. Earlier, Abrahams and his collaborators carried out extensive simulation
studies of the rare gas atoms physisorbed on graphite [9]. These simulation studies
were on LJ systems. In this chapter, we use Molecular Dynamics (MD) simulations
to study the role of incommensurate corrugation potential for a 2D system of atoms

interacting via a repulsive Yukawa potential, a realistic representation of intercala-

 

89

tion compounds where the intercalants are charged objects. It should be noted here
that screened Coulomb systems, such as the one addressed here, take significantly
more computation time (due to the exponential form of the interaction) compared
to the LJ and other (algebraic) “short range” systems. Although we expect some
general characteristics of Lennard-Jones and repulsive Yukawa systems to be similar,
it is of interest to compare the detailed structure and the nature of the transition

in these two potential systems [18].

The arrangement of the chapter is as follows. In Sec. 5.2 we introduce three
different systems. The temperature dependence of different physical quantities are
studied for these systems in Sec. 5.3. Finally, in Sec. 5.4 we summarize the important
results of the present work and discuss its implications vis-a—vis earlier work on LJ

and related systems.

5.2 Physical Systems

Details of potential of the system and MD simulations were discussed in chapter 4.
Here, we introduce three systems which will be studied for 2D melting. MD simula-
tions have been performed with 216 Rb ions distributed to ensure a planar density
of 12 Carbons/Rb (corresponding to the stage-2 Rng4 system), i.e., 0.0318 Rb
ions/A2. We call this system as system II. Starting from this system, we construct
two other systems. System I (2K = 0 in Eq. (4.1)) has the same density as system
11 but without any corrugation potential. By Comparing results of the systems I
and II, we can study the effects of corrugation on various physical properties. The
ground state structure of the system I is an equilateral triangular structure. In
addition to these two systems, we have also studied a 301-particle system (system

III) with a slightly smaller Rb density (0.031] Rb ions/A2) which corresponds to

a

90

the stage—2 Rngu-z stoichiometry. The ground state of the system III turns out to
be a perfectly ordered periodic domain wall (PDW) structure [19] where hexagonal
placquettes consisting of 7 Rb ions (nano domains) are periodically arranged, sepa-
rated by a periodic array of domain walls (see Figure 4.3(b)). The ground state of
III is commensurate whereas that of II is incommensurate. Note however that even
in the commensurate case, only the atoms at the center of the J7 x \/7 domains
sit at the minima of the single particle potential and the other atoms are displaced
from these minima, but very slightly [20]. Various types of PDW structures were
extensively discussed in chapter 4. The ground state of RbC'u (System II) can be
described by the above PDW structure with additional (about 2.2% more) Rb ions
occupying interstitial sites (see Figure 4.3(a)). We will refer to these additional
Rb atoms as atomic (distinct from topological) defects. The effect of these atomic
defects on physical properties, particularly low energy excitations will be discussed
in chapter 7. The reason for choosing these two systems with close densities is to
see whether the melting process and related thermodynamic quantities are sensitive

to the atomic defects.

5.3 Thermodynamic Properties

We discuss the melting of the systems I] and III and compare their melting pro-
cess with that of a corrugation free system (system I), the usual triangular solid.
Temperature dependence of energy, bond orientational order parameter, its angu-
lar susceptibility, topological defects, and translational diffusion constant will be

studied to monitor melting of these three systems.

91

5.3.1 Energy vs. Temperature

In Figure 5.1, we give the temperature dependence of the potential energy (total
energy-IcBT) per particle as a function of T both for the heating and the cooling
runs. In the absence of the corrugation potential (2K = 0), there is a rapid change
in this energy in the temperature range 195K < T < 200K. The transition is most
likely first order as seen in previous MD simulation of other physical systems (LJ
systems [9, 21], 1/r" systems [22] and e‘r" / r system [23], Weeks-Chandler-Anderson
(WCA) system [24] which is a truncated version of the LJ system). The absence of
hysteresis might be due to either the small size of our systems or weak first order
nature of the transition. We note here that the earlier MD simulation studies of
a repulsive screened Coulomb were carried out for smaller (100 particles) systems
in the absence of corrugation and showed hysteresis effects [23]. The details of the
equilibration time at different temperatures in this study are not available but it
is possible that the observation of hysteresis might be the result of nonequilibrium

effects.

If we assume that the system I shows a first order transition and take the slopes
of the energy vs T above 210K and below 190K, then we find the change in entropy
associated with this first order transition to be AS/kg = 0.23 which agrees rather
well with the earlier simulation results on other systems. As seen in Figure 5.1, the
effect of corrugation is drastic; in contrast to the corrugation-free case (system I), the
potential energy changes smoothly for both systems II and III in the temperature
range 190K < T < 230K, although there is some indication of a slope change near
220K. To make sure whether there is actually a solid-liquid phase transition in
systems II and III, we must calculate other physical quantities which are sensitive

to this transition. In addition, one must increase the strength of the corrugation

92

 

T'lIIITIITI'IIIVVIIII'

—

if“
1- *#*xx
- I

p
‘1

p
G
T

.‘f
21
at
‘9.

TE — KE (1n umts of 300 kg)

in...
a

 

 

lllllLlllllllllllllllll-i
100 150 200 250

T

 

Figure 5.1: Potential energy vs. temperature with and without corrugation poten-
tial. Energy is measured in units of 300193 and the temperature in units of Kelvin.
System I corresponds to RbCu without corrugation (2K =0), system II corresponds
to Rng4 with corrugation (2K=0.9), and system III corresponds to Rng4,57 with
corrugation (2K=0.9). 2K is a measure of the strength of the corrugation potential
(see Eq. (4.1)). For all these three systems both cooling and heating results are
shown.

93

systematically to see how the corrugation alters the nature of the melting transition.
5.3.2 Bond Orientational Order Parameter

Since it is difficult to investigate the melting transition by calculating the posi-
tional correlation function for small systems, one usually monitors the temperature
dependence of the bond orientational order parameter (BOOP) and its associated
susceptibility. Strandburg, Zollweg and Chester [25] introduced the BOOP by defin-
ing a quantity the as

88°60”

 

¢6 =71v- s (5.1)

 

 

(j "I
where the sum on I is over all the N particles, and the sum on j is over the nearest
neighbours of a particular particle I, n; being the number of such nearest neighbours

defined by the standard Voronoi construction. The angle 0;,- is the orientation of

the bond joining the particles 1 and 1' measured with respect to a given fixed axis.

In Figure 5.2 we plot #2.; as a function of T for systems I and II. In the corrugation-
free case there is a rapid decrease in the order parameter between 190K and 200K
whereas in the presence of corrugation the change in the order parameter appears to
be less abrupt and the transition temperature is somewhere between 210K and 230K.
In addition, the order parameter in the solid phase is larger in the corrugation-free
case. This can be easily understood by recalling the ground state structure of these
systems. The system without corrugation has a perfect triangular lattice structure
at T=0K for which ¢6=1 whereas in the system with corrugation the domain walls
reduce the strength of this order parameter. The nonzero value of the order pa-
rameter in the liquid phase of system I is due to finite-size and periodic boundary
condition because 1/16 should vanish in the thermodynamic limit. In the presence

of the corrugation potential of hexagonal symmetry 31).; should be nonzero even in

EPA.

94

 

1.0 IIIITIITMUIUI'TWFI'Y‘IIT'

1

(L8
(L6
(L4
0.2 X. 0 : Heating K

+. D : Cooling <> 0

0.0 Lillllllllljlllljlelll I]
100 150 200 250 300

T

lllllllllllllJllllMllLJ

“We
TrYUlU'TTIjjjiliii'Ii‘lfT

 

 

 

Figure 5.2: The bond orientational order parameter 1126 as a function of temperature
for systems I and II (see Figure 5.1 captions for explanations). Both heating and
cooling run results are shown.

95

the thermodynamic limit. We however note that in our finite system simulations
(with periodic boundary conditions), for temperature above 240K the residual or-
der parameter is independent of the strength of the corrugation suggesting that the

influence of finite size effects in the calculation of 1123 is large.
5.3.3 Angular Susceptibility

To further probe the nature of the transition and to obtain the transition tem-
perature more accurately we have calculated the susceptibility associated with the

BOOP. We define the angular susceptibility x5 as
M = N [< [¢6[2 > -- < label >’]/kBT. (5.2)

In Figure 5.3, we give the temperature dependence of x5 for the two systems dis-
cussed in the previous paragraph. In the corrugation-free system, the susceptibility
increases rapidly with decreasing T and shows a dramatic (almost discontinuous)
dr0p between 200K and 195K. This is the same temperature range where the order
parameter shows a rapid change (Figure 5.2) thereby suggesting that the transition
is most likely first order. From our finite size simulation we cannot conclusively
state that the transition is indeed first order. In the absence of an observable hys-
teresis we cannot rule out a sharp but continuous transition. Glasser and Clark [24]
in their extensive MD simulation studies of a truncated LJ system found a nearly
discontinuous drop in x6 as a function of increasing density (temperature being
held constant) indicating a first order transition. They however did not discuss the

hysteresis issue.

In the presence of corrugation the behaviour of the susceptibility is much more
intriguing. For the Rng, system (system II), X6 is relatively small and is a smooth

function of temperature showing a very broad peak near 230-240K. There is not

96

 

IIIIIIUTITIUIIIII'I'IT‘III

P I

 

 

 

 

. + .
6— + -—
C x i

+
_ I ++ .

m +
5‘ 4” if "
g . .

co ’$<
* p -
>3 - 2x .
b x -l
l— _

2 i +

_ II» .
“$32939“ *
o LII 0L11210191L|llll|11
250 300

T00

Figure 5.3: The angular susceptibility x6 (see Eq. 5.2 for the definition) vs. tem-
perature for systems I and II (see Figure 5.1 captions for explanation). Heating and
cooling run results are shown.

97
much action near the temperature where the order parameter shows a rapid drop
at 200K< T (220K (indicated by an arrow). Clearly the melting transition of this
domain wall solid is very different from that of a simple triangular solid (system
I). In order to pin down the precise nature of the transition it will be necessary to
carry out extensive Monte Carlo simulations in larger systems covering the transition

region and also the region where the susceptibility shows a peak.

One plausible physical argument we can give to explain the temperature de—
pendence of 11’s and x3 in the domain wall solid case is the following [26]. As we
increase the temperature, the region near the domain walls starts to soften or ‘melt’
(~ 200K); this leads to a rapid decrease in the long range bond orientational or-
der. However above this transition there is still sufficient short range orientational
correlation inside each triangular domain and the peak in x6 can be thought of
as a second ‘melting’ of these highly correlated small triangular lattice clusters of
J7 x J7 domains, the transition region being broadened by finite size effects. The
fact that this broad peak is seen at temperature ~ 240K which is higher than 200K
where the corrugation-free J6 x 1,56- triangular solid melts can be understood by
the following simple argument. Ordinarily a J7 x J7 triangular solid in the absence
of corrugation will melt at a temperature lower than 200K because of lower energy
compared to the J6 x J6 case whose melting temperature is about 200K. However
in the presence of a strong corrugation potential the melting temperature will be
increased from this value and it is possible that the melting temperature is near
240K where one sees the broad peak in x6. In fact, Vives and Lindgard [16] in their
MC simulation study of a commensurate LJ system found an increase in the melting

transition temperature with an increase in the strength of the corrugation.

[0"

98

Since the ground state of My can be described by a perfectly ordered PDW
solid with a few additional Rb atoms (defects) and since it is reasonable to expect
that the order parameter, the susceptibility and the nature of the transition might
be affected by these defects, we have studied the temperature dependence of tbs and
x6 for the stage-2 Rng4,57 (system III) whose ground state structure is shown in
Figure 4.3(b). In fact, one can study the melting properties of a whole class of PDW
solids where the widths of the domains and domain walls can be systematically varied
by simply changing the planar Rb concentration [19]. This can in principle tell us
how the nature of the melting transition changes as we introduce domain walls into
the system. For Rng4,57 where the J7 x J7 triangular domains contain 7 atoms
separated by 2 x 3 x J7 walls, the order parameter at low temperatures is slightly
larger than RbC'u [see Figure 5.4(a)] consistent with the presence of additional
defect atoms in the latter system. For the former we see a phase transition (rapid
decrease of 1,126) near T=230K [Figure 5.4(a)] with a peak in the susceptibility x3
[see Figure 5.4(b)] at a higher temperature (T ~ 240K). This is qualitatively similar
to that shown by the system II. Thus we conclude that the peak in orientational
susceptibility above the melting transition temperature appears to be a general
property of the melting of domain wall solids. The removal of defects makes the
susceptibility peak in the system III slightly more pronounced compared to the
system II. But in both cases the peak intensity is about 3-5 times smaller than the

corrugation free case.

5.3.4 Topological Defects

To further explore the effect of corrugation on melting we have studied the tem-
perature dependence of the concentration of local topological defects (LTDs) which

are bound pairs and quartets of 5- and 7-fold coordinated atoms both in the pres-

99

 

.(..)I...iliiiiTwiiilii.iliiq
a
0.8
III
0.6 - r T I ’ f‘:‘:a"dl’f"}:7.;l7=:%
; -
0.4 II

C
N
gjllllllllllllelLl

 

 

0.0 ,
2 o (b) -E
a: III 3
1 5 W ~—
8 pt 5
co >63 % _-
3'; 1.0 cg .
>‘<’ si . Wi
0.5 x a“.
I '. j
000 $173k“. I" H.151? .' .210! l l l l l IMLJ L l "i

100 250 300

T00

Figure 5.4: Temperature dependence of (a) bond orientational order parameter tbs
for systems II and III and (b) corresponding angular susceptibility Xs- See Figure 5.1
captions for the description of different systems. Both heating and cooling runs are
shown.

100

ence and the absence of corrugation (see Figure 5.5). We use the standard Voronoi
polyhedra construction to obtain the number of 5, 6, and 7 fold coordinated atoms.
It appears that the effect of corrugation on the density of these LTDs is marginal
excepting in the transition region. In the transition region, the temperature depen-
dence of the LTD density pp is smoother for the corrugation case (see Figure 5.6)
indicating a continuous transition. One interesting observation is that in the solid
phase (T < 195K) and in the liquid phase (T > 250K) pp is weakly dependent on
the corrugation strength (2K). The corrugation potential has a strong effect on the
defect density only in the transition region. When 2K =0, we see that the rate of
increase of pp with temperature, dpp/dT, is maximum at T=200K which is also
the maximum of x,;. In contrast, for 2K = 0.9 i.e. for the system II, dpp/dT is
maximum near T=220K whereas x6 peaks at T ~ 240K, again indicating a basic
difference between the melting characteristics of a triangular solid in the absence of

corrugation and those of a domain wall solid in the presence of corrugation.
5.3.5 Translational Difi'usion Constants

Finally, the translational diffusion constant D using the relation (1") = 4Dt has
been calculated to monitor the solid-liquid transition and the effect of the corru-
gation on the diffusion rate. In contrast to its effect on the density of topological
defects, corrugation profoundly affects the diffusion rate (see Figure 5.7). The nearly
first order nature of the transition in the corrugation free case is seen as a dramatic
increase in D at about 200K. On the other hand, in the presence of corrugation D
starts to increase rather slowly near the transition temperature (T=220K), seems
to flatten slightly near 240K where there is a peak in the orientational susceptibil-
ity, and finally increases linearly with T. The drastic reduction in D in the liquid

phase in the presence of corrugation can be understood in terms of an additional

101

 

 

 

 

 

Figure 5.5: Topological defects (5- and 7-fold coordinated atoms) (a) in a solid and
(b) a liquid state obtained by using Voronoi polyhedra construction.

102

 

 

 

100 . ' I I I] I I I I I I I T I l V—I y 1 _,
:9 —:
b 0 ‘

*3 3 813%“ 1
.2 60 " I §8 -
9’ F 4* 4
“a . , Bo .
4O — _.

‘6 I. % D H "
no : X0 .
5 ' a 3
20 — ._

z i a? :
1' .

ON 1 l 1’ 1 ’ ’ ' I ‘ ' 1 1 l I 1 1 1

 

100 150 200 250 300
T

Figure 5.6: The number of disclinations (topological defects, 5- and 7-fold coor-
dinated atoms) vs. temperature. Observe the sharp increase in the number of
disclinations near 200K for system I (without corrugation) and smooth increase
near 220K for the system II (with corrugation). Different symbols refer to heating
and cooling runs.

I»—

103

 

4 I I I III I I I l I I rTW
,. x. .
3— —
fa . a
a - + + .
N\ . X 4
E ' f i
O 2— + _
"1’ - 3* -
C - I + -
H
v 1- .-
Q - X .
1— —

 

 

 

 

Figure 5.7: Diffusion constant D vs. temperature for system I (without corrugation)
and system II (with corrugation). Both heating and cooling runs are shown.

104

activation process involving the corrugation potential. Using the parameters of our
potential we estimate the energy barrier associated with the corrugation potential
to be 1080K. The experimental value of the activation energy obtained from the dif-
fusion measurements is ~ 0.063 eV (756K) which is about a factor 2/ 3 of the single
particle activation barrier. The net activation energy is of course a combination of
the single particle and interaction potential. More detailed temperature dependent
study of the diffusion rate D in the liquid phase is required to quantitatively deter-
mine the true activation energy. Further details of the diffusion and the nature of

the dynamics will be discussed in chapter 6.

5.4 Summary and Discussion

In summary, we have carried out extensive MD simulations of the melting and
freezing transitions of a 2D repulsive Yukawa system in the presence of (strong)
incommensurate corrugation potential. This model describes rather well the physical
properties of stage n (n 2 2) graphite intercalation compounds and similar systems
where the intercalants are charged objects interacting with a screened Coulomb
potential. The parameters we have chosen are appropriate for stage-2 Rb GIC.
Before comparing our simulation results with experiments in stage-2 Rb GIC we

summarize some of our general observations.

The triangular lattice structure of the corrugation-free solid changes to a periodic
domain wall structure in the presence of a (strong) incommensurate corrugation
potential. The size of the domains and domain walls depend on the strength of the
corrugation potential and the intercalant density. Whereas the former (triangular
solid) undergoes a sharp (perhaps discontinuous) melting/freezing transition, the

latter (PDW solid) shows a smooth transition. The transition temperature increases

105

in the presence of corrugation. Such an observation was also made by Vives and
Lindgard [16] in their Monte Carlo study of the melting of J3 x J3 commensurate LJ
systems. In the absence of corrugation, we find that the rapid decrease in the bond
orientational order parameter p, and the peak in the corresponding susceptibility
x3 occur at the same temperature. In contrast, in the presence of the corrugation
potential, )0; does not show any structure (discontinuity or peak) at the temperature
where 4).; shows a rapid decrease, but has a rather broad weak peak several degrees
above this temperature. We have given a physical explanation of this behaviour in

terms of the melting of a domain wall solid.

The temperature dependence of x6 for the corrugation-free system shows a sharp
A-like peak. The ratio of the half width AT to the transition temperature T i.e.
AT/T=0.l. Recently Glasser and Clark [24] have carried out extensive MD sim-
ulations on a truncated LJ system (WCA system) and have studied the melting
transition at a constant T but as a function of density (p). They also find a sharp
change (actually a discontinuity) in x5 near the melting transition indicating a dis-
continuous transition, very much like what we have seen as a function of temperature
at constant density. The half width of their peak in x3 vs p, Ap/pc=0.22. In con-
trast to our MD results, the MC results (for LJ systems) of Vives and Lindgard [16]
in larger systems (N=2700) show an extremely broad peak in x3 and appreciable
hysteresis, most likely due to nonequilibrium effects. Although our system size is
small, our results are for well equilibrated systems as indicated by near agreement
of the heating and cooling runs. Because of the above difficulties we are unable to
compare the sharpness of the susceptibility peak (as a function of temperature) in

the corrugation free limit for repulsive Yukawa and LJ systems.

Finally we would like to compare our simulation results with available experi-

106

ments. It is not possible to directly measure the bond orientational order parameter
and the associated susceptibility. The former can be indirectly measured from the
anisotropy of the X-ray (or neutron) diffraction peaks in the “solid phase”. Zabel
and coworkers [27, 28] measured the intensity of the superlattice peak as a function
of temperature using neutron diffraction. They find a continuous transition from
the solid to the liquid phase at about 165K in stage-2 RbCu. Our simulation re-
sults agree with this except that our transition temperature is about 220K. This is
due to (1) finite size effect and (2) limitations in the potential model that we have
used. In fact, Moss and coworkers [29, 30] have used a different 2-body potential
(but the same l-body potential as used here) and find the transition temperature
to be somewhere near 160K. However they did not carry out a careful T-dependent
study of different physical quantities as we have done here. By comparing the two
different 2-body potentials in the region of physical interest we find that the Moss
et al.’s potential is slightly weaker than ours. This might explain the difference in

the transition temperature discussed above.

Bibliography

[1] J. M. Kosterlitz and D. J. T houless, Progress in Low Temperature Physics, Vol.
VII-B, ed. by D. F. Brewer (North-Holland, Amsterdam), p.373.

[2] D. R. Nelson and B. I. Halperin, Phys. Rev. B 19, 2457 (1979).
[3] A. P. Young, Phys. Rev. B 19, 1855 (1979).
[4] K. J. Strandburg, Rev. Mod. Phys.60, 161 (1988).

[5] R. J. Birgenau and P. M. Horn, Science 232, 329 (1986), N. Grieser et. al., Phys.
Rev. Lett. 59,1706 (1987).

[6] Q.M. Zhang and J. Z. Larese, Phys. Rev. B 43, 938 (1991).
[7] s. T. Chui, Phys. Rev. Lett. 48,933 (1982).

[8] W. A. Steele (Private communication).

[9] F. F. Abraham, Adv. Phys. 35, 1 (1986); F. F. Abraham, Phys. Rev. B 28,7338
(1983); ibid B 29,2824 (1984).

[10] AD. Novaco and J.P. McTague, Phys. Rev. Lett. 38, 1286 (1977).

[11] L. K. Runnels, Phase Transitions and Critical Phenomena, ed. by C. Domb and
M. S. Green, (Academic Press, New York, 1972), Vol 2, p.305.

107

108

[12] M. J. Colella and R. M. Suter, Phys. Rev. B 34, 2052 (1986).

[13] K. L. D’Amico, J. Bohr, D. E. Moncton, and D. Gibbs, Phys. Rev. B 41, 4368
(1990).

[14] J. Cui and s. C. Fain, Jr., Phys. Rev. B 39, 8628 (1989).

[15] S. C. Moss and R. Moret, Graphite Intercalation Compounds, Vol 1, Structure
and Dynamics, ed. by H. Zabel and S. A. Solin, (Springer Series on Topics in

Current Physics, Springer (1990).
[16] E. Vives and Per-Anker Lindgard, Phys. Rev. B 44, 1318 (1991).
[17] E. Vives and Per-Anker Lindgard, (private communication).

[18] Recently structural phase transition (from BCC to FCC) and melting transition
for a 3dimensional repulsive Yukawa system has been studied extensively by M.
0. Robbins, K. Kremer, and G. S. Grest, Jour. Chem. Phys. 88, 3286 (1988).

[19] Hyangsuk Seong, S. Sen, T. Cagin, and S. D. Mahanti, Phys. Rev. B 45, 8841
(1992).

[20] The structure of the PDW solid is quite similar to the one proposed by J. Cui
and S. C. Fain, Jr. (see Ref. [14]) for the A-phase of Dg/graphite where J3 x J3
domains and J3><1 walls are periodically arranged. One should also note that all

the D2 molecules do not sit at the hexagon centers as in the PDW model (see

Ref. [17]).
[21] C. Udink and J. van der Elsken, Phys. Rev. B 35, 279 (1987).

[22] P. Vashishta and R. K. Kalia, Melting, Localization, and Chaos, ed. by. R. K.
Kalia and P. Vashishta, North Holland (New York) 1982, p.43.

7’2.-

109

[23] H. Chen, P. Dutta, D. E. Ellis, and R. Kalia, Jour. Chem. Phys. 85(4), 2232
(1986).

[24] M. Glasser , Ph. D thesis, University of Colorado (1991).

[25] K. J. Strandburg, J. A. Zollweg. and G. V. Chester, Phys. Rev. B 30, 2755
(1984).

[26] There is an alternate way of looking at the melting transition (suggested to us
by Prof. S. Jayaprakash of Ohio State University) i.e. in terms of condensation
of defects which are in this case domain walls. One can visualize preexisting
elemental J7 x \f7— triangular placquettes in the liquid phase above the solid-
liquid transition. The region between these placquettes (domain walls) can be
thought of as defects which then order below a critical temperature. For the
system size we have studied, the number of these defects is small (by about a
factor of 7) compared to the number of atoms and therefore one sees a rather
broad transition region. To test this interesting idea one has to study the melting
of the PDW solid for larger systems and see if the transition region sharpens with

the size of the system.

[27] H. Zabel, A. Magerl, J. J. Rush, and M. E. Misenheimer, Phys. Rev. B 40,
7616 (1989).

[28] J. D. Fan, George Reiter, and S. C. Moss, Phys. Rev. Lett. 64, 188 (1990).
[29] S. C. Moss, et al., Phys. Rev. Lett. 57, 3191 (1986).

[30] J. D. Fan, Omar A. Karim, G. Reiter, and S. C. Moss, Phys. Rev. B 39, 6111
(1989).

Chapter 6

Dynamics of 2D Repulsive
Screened Coulomb Fluids on a
Corrugated Surface

6.1 Introduction

Excepting for a short discussion on the diffusion constant, we have been until now
concerned with equilibrium (static) properties of 2D systems. Next two chapters
will focus on the question of dynamics. Dynamics of homogeneous fluids in three
dimension is well understood [1]-[3]. In recent years there has been considerable
interest in understanding the dynamic (and thermodynamic) properties of fluids
subjected to different types of external constraints [4]-[11]. Examples of such sys-
tems are: two dimensional fluids in the presence of external corrugation potential as
in intercalation compounds [6]-[9], premelted surface fluid phase [10], and fluids in
porous media [11]. One expects the particles in the fluid phase to exhibit different
types of dynamic behavior in different length and time scales due to the inhomo-
geneous nature of these external constraints. An extreme example is the solid-like
dead layer formed near the constraining walls when a fluid flows through a meso-

porous medium. Here one expects to see the coexistence of solid-like and fluid-like

110

111
dynamics.

Dynamic properties of physical systems can be studied by time correlation func-
tions and their Fourier transforms which can be directly obtained from neutron or
light scattering experiments. In neutron scattering experiments thermal neutrons
are directed onto the fluid with an incident wave vector k; and energy E.- and the
scattered neutrons are measured as a function of the wave vector k, and energy
By. One measures the scattering cross section (a) as a function of q = k, — k; and
Iw = E f - E5. The cross section directly relates to the dynamic structure factor
S(q, v) [1, 7]. Even if the fluid is in a state of macroscopic equilibrium, spontaneous
microscopic fluctuations in the local density occur in the system due to thermal ex-
citations, and these fluctuations give rise to the above dynamic scattering. Usually
in thermal neutron scattering [12] from fluids, density fluctuations in space and time

lead to characteristic q and u dependence of S (q, V).

In this chapter we address the dynamics of a specific class of layered systems,
namely 2D fluids subjected to a periodic corrugation potential whose strength can
be tuned. Physical realizations of such systems whose dynamics has been extensively
probed experimentally are graphite intercalation compounds (GICs) [6]. Zabel and
coworkers [7, 8] have carried out neutron scattering experiments in alkali metal
GICs for which the substrate corrugation potential produced by the graphite layer
on the alkali intercalants is strong. However, to study the evolution of the dynamics
from a homogeneous fluid to a lattice-fluid, we change the strength of the periodic
corrugation potential from 2K = 0 to 2K = 0.9 (see Eq. (4.1) of chapter 4 for
definition of 2K) through intermediate values, 0.15, 0.36, 0.55, and 0.75, the value
of 0.9 appropriate for stage-2 Rb GIC. See chapter 4 for the the details of potentials

and Molecular Dynamics (MD) simulations. Here we study larger systems with 864

112

particles and concentrate on the properties of the liquid state obtained by averaging
over the directions of the scattering wave vector in the calculation of correlation

functions.

To understand the dynamics of fluid systems in detail, we study different cor-
relation functions such as van Hove correlation function, velocity auto-correlation
function, intermediate scattering function, and dynamic structure factor. These will
be defined in the following sections and we will discuss how these functions are af-
fected by changing the strength of the corrugation potential (i.e. in going from a

simple homogeneous liquid to a GIC system).

6.2 Time Correlation Functions in Real Space

An equilibrium density-density time correlation function was first introduced by van
Hove [13] to discuss the density fluctuation in a liquid. It is given by
1 N N
G(r, 1) = N<§j§§ 5(1' + l“5(0) - 17(0)), (6-1)
where r,-(t) is the position of the i-th particle at time t and the angular brackets
imply time average. Physically, G (r, t)dr is proportional to the probability of finding
a particle i in a region dr around a point r at a given timet that there was a particle

j at the origin at time t = 0. The function G(r, t) can be separated into two terms,

usually called the “self” (8) and “distinct” ((1) parts, i.e.
G(r,t) = G,(r,t) + Gd(r, t), (6.2)

where

N
0.9.9 = 71,430+ r40) - mm». (63)

5:1

113

N .
Gd(r, t) = fig: 6(r + r,-(0) — r.-(t))>- (6.4)

Rom the above equation we find that G.(r,0) = 6(r) and Gd(r,0) = pg(r), where
g(r) is a pair distribution function which gives the probability of finding a pair of
atoms a distance r apart. For isotropic fluids, both 6'. and 04 are functions of the

scalar quantity r.
6.2.1 G,(r, t) for Typical Systems

Before discussing the space and time dependence of G.(r, t) for a corrugated fluid,
let us review what is known about this quantity for simple physical systems. For an
ideal gas, a diffusing atom, and a harmonic oscillator, it is found that the self part
of the function is a Gaussian in r [1, 14]. Furthermore, it is rigorously true for any
dynamical system as t -+ 0, because the atoms behave as if they are free. Therefore,
one usually assumes that the spatial dependence of G.(r,t) is approximated by a
Gaussian [12],

 

1 "2 2
__ -r /2a(t)

where d is the dimension of the system. The mean square displacement gives a

physical meaning to a(t), i.e.,

((r(t) — r(0))2) = / drr’G,(r, t) = a(t)d, (6.6)

namely a(t)d is the mean square displacement after a time t for d-dimensional
system [see Appendix C for the proof of Eq. (6.6)]. There are some cases for which

a(t) is known exactly. For free particles using r(t) — r(0) = tp/m, one finds that

a(t) -_- “ET“, (6.7)

m

 

114

where T and p are respectively the temperature and the momentum of the free

particle. For a diffusing Brownian particle of mass m
Dm
(1U) — 2D(t - $71), (6.8)

which is obtained easily from the equation of Brownian motion (i.e. Langevin equa-
tion), D being the diffusion constant [14]. For harmonic phonons in solids, one

obtains

2kBT

-coswt

m]... de(w )— (6.9)

 

0(1)-—

where Z (w) is the normalized density of phonon states. Eq. (6.9) recovers the free
particle result when we expand cos wt in powers of small t and keep the leading
term. However, at long times this does not increase continuously but approaches a
finite limit. In fact, this particular feature distinguishes a solid from a liquid and a

gas.

We know explicit forms of a(t) for short and long times for a simple liquid.
However, for intermediate times considerable deviation from a Gaussian behaviour
occurs. Except for small times, the Gaussian form also does not hold for a particle
diffusing in a solid where the particle stays mainly at discrete positions determined
by the lattice. Egelstaff and Schofield [14] constructed a simple a(t) for intermediate
time t which recovers the correct short and long time behavior for a liquid. The

expression for a(t) is given by
a(t) = 2D[(t2 + c2)”2 — c], (6.10)

where c is a parameter characteristic of the model. An easy way to study the

dynamics for intermediate times is through MD simulations.

115

6.2.2 Self-Diffusion Constant

A macroscopic quantity which characterizes the diffusion processes in a one-
component fluid is the self-diffusion constant D. By looking at the mean square
displacements of particles, one can study the diffusion process. Figure 6.1 shows
trajectories of a tagged particle for various strengths of the corrugation potential
but for a fixed temperature. The figure clearly shows that the motion of a particle
is confined as the strength of the corrugation potential increases. In a given time
interval, particles of the 2K = 0 system can diffuse much farther than those of the

2K = 0.9 system. For a 2D system, the diffusion constant D is defined as,

 

2
D - lim (Ar) .

— t-ooo 4t

(6.11)

To see how D approaches its t —v 00 limit, we have calculated a time dependent

diffusion coeflicient D(t)

D(t) = 94:12. (6.12)

Figure 6.2 shows the time dependence of the D(t) for various strengths of corrugation
potential near 250K. Whereas for the 2K = 0 and 0.15 systems D(t) increases very
slowly, for the rest of the systems D(t) decreases even more slowly as time goes on.
These different trends might be understood by looking at the barrier height A (in
units of kg) of the corrugated potential. For T near 250K, the first two systems are
at temperature higher than A (0 for 2K = 0 and 180K for 2K = 0.15), whereas
A’s for the other systems are larger than 250K (see Table 6.1). This t-dependence
of D(t) might tell us something about the anomalous behaviour of diffusion in 2D
Systems resulting from the long time tail of the velocity auto-correlation function

[15]. To clear up this question, we need to study ( 1) larger size of the systems so

116

 

TIUW UIIWT‘W'IIUjltlU'IUIIVI

L

2K=0 2K=0. 15 2K=0.36

We

2K=0.55 2K=0.75 2K=0.9

51’ as e

_16lllllJllllllllljJJLlllllllll

--15 -10 -5 O 5 10

10

O a
IIIIITirYIYIUUIYYIIIIVfiIYII

lJlllllllllllllllllJlJllllJ

 

 

 

Figure 6.1: Trajectories of tagged particles during about 100ps for various 2K at a
temperature near 250K.

117

 

0.6 _ ............................................. _

 

...
Q.-
o.-.
0.0-.0-0000-0...-coon-oooo--IOOO-....

...~
...-
..'-..
0--..o-ooo-ooo-coo-ooo_...-...-..

 

 

Figure 6.2: Time dependent diffusion coeflicient at about 250K. From top to bottom
curve, 2K = 0, 0.15, 0.36, 0.55, 0.75, and 0.9.

 

118

Table 6.1: Barrier height (A) from the various corrugation potentials, the temper-
ature (T) of corresponding systems, and the diffusion constant (D) obtained from
mean square displacement (and D obtained from velocity auto-correlation function)

for different 2K.

 

 

 

2K A T (Kelvin) D(10’4cm2/ sec)
0 o 254 0.633 (0.616)
0.15 180 253 0.479 (0.479)
0.36 432 251 0.253 (0.254)
0.55 660 254 0.141 (0.141)
0.75 900 256 0.062 (0.070)
0.9 1080 248 0.033 (0.036)

 

 

 

 

 

 

that our (r2) calculations are not affected by periodic boundary conditions and (2)
carry out longer time runs to see this ‘long’ time behavior. However, in our present

MD simulation studies, we do not approach this long time behavior.

Diflusion constant also can be obtained from the time integral of the velocity
auto-correlation function (VAF). The VAF is not directly measured in scattering
experiments and the most direct method of determining the VAF is actually MD
simulations. It reveals molecular motions on a microscopic time scale where dynam-
ical details occurring over intervals of ~ 10‘13sec in real time can be resolved. Since
the interval is comparable to the molecular collision time, one is effectively probing

features of dynamical behavior most sensitive to the details of the interactions.

The VAF is related to the diffusion constant as follows,

(Ar)2 = 2th (6.13)

= fo'vu'w . jtv(t”)dt” (6-14)

0

119

= 2t [0' ds(v(s) . v(0)) — 2 [0' sds(v(s) . v(0)), (6.15)

where d in Eq. (6.13) is a dimension of the system. In the limit of large t, we get

Dt=:11-l_1.m/0ds(v(s)- v(0)). (6°16)

The 4-th column of Table 6.1 shows diffusion constants for various strengths of the
corrugation potential. The first numbers obtained from the mean square displace-
ments agree very well with the second ones (given inside the parenthesis) obtained
from the integration of the VAF. To carry out the numerical integration accurately
in Eq. (6.16), each point is joined smoothly and the interpolated values are used in

the integration.

In addition to the diffusion constant obtained from the time integral of the
VAF, the latter’s time dependence conveys a large amount of information about the
dynamics of a particle in the real time domain. To understand this t-dependence,
let us study the VAF in more detail. The normalized VAF is defined by

(E.- V10) ° W(O»

Z“) = (2:.- v.-(0) -v.-(0)>'

 

(6.17)

Figure 6.3 shows the VAF for various values of 2K near 250K. As in an ordinary
liquid, Z (t) shows an oscillatory t-dependence with decreasing amplitude resulting
from the diffusive dynamics. The effect of corrugation at least up to t ~ lps is
to decrease the area under the positive part and increase that under the negative
part. This makes a major contribution to the overall decrease in D with increasing
corrugation strength. The characteristic time dependence of Z (t) can be described
broadly with three regimes, short time, intermediate time, and the long time be-

havior.

 

120

 

TWIIIUTIIIUIIIIIIIIITIV

 

1.0 -'
' 250K .
q
l- ..
I- d
0.5 l '4
A " 1 "
no.) a
v " o . "
N . k
'- I :’ ‘0' J.”- "‘
I i ' °- . ‘ ,1
000 _ .l ‘3 .”:. o..‘\o‘fi.\ T.“
f on \fi" '1
.. :r.‘ ' .
.\_,
I- 2‘ '1
P .‘.. all
-0.5 F ' '-

 

 

plllmllillllllllllllllLll.

0 1 2 3 4 5
t (138) A

 

Figure 6.3: Velocity auto-correlation functions for various 2K near 250K. From the
solid (bottom) curve, 2K = 0.9, 0.75, 0.55, 0.36, 0.15, and 0.

121

Short-time behavior of the VAF can be studied by Taylor series expansion, i.e.
v,(t) = v,-(0) + v,(0)t + £940)? + - .. (6.18)

Multiplying by v,(0) and taking a canonical ensemble average we obtain

<v.-(t)-v.-(o» = <v.-’>+§<v.-v.->t’+-~ (6.19)
= (v?)—-;-(13.-2)t2+~- (6.20)
= <v?>(1-§93t’+-~). (6.21)

where {23 = (|F|2)/2mchT. Time reversal symmetry makes the odd terms in t
vanish. The short time VAF is related to the mean square acceleration (or force)
which is one of the basic properties of the VAF and provides information about
the interaction between particles, single particle potential, and certain structural
properties of the fluid. Clearly increasing the strength of corrugation potential

increases (13 resulting in a faster decrease in Z (t) starting from Z (0) = 1.

For homogeneous liquids, it is known that the intermediate time behavior of Z (t)
is quite sensitive to the temperature and the density. At low density, VAF usually
decays monotonically which implies the absence of many body correlation effects. As
the density increases, a minimum in the VAF occurs indicating a resonant behavior
or the presence of memory effects (which produce correlations between random
force and velocity) in the system. Z (t) for hard spheres decreases monotonically as
a function of t, with a weak but well-defined tail up to about thirty mean collision
times. At higher density, VAF shows a strong oscillation qualitatively reminiscent of
a damped oscillator and takes on negative values. This behavior can arise from the
tagged particle being confined in a cage formed by its immediate neighbors. Negative
Z (t) means that at timet a particle moves in opposite direction to that at timet = 0.

Therefore, Z (t) describes a “back scattering” effect, a result of collisions between

 

122

nearest neighbors which leads to a reversal of velocity into a relatively narrow range
of angles. There is a long time tail controlled by hydrodynamic conservation laws,

but we are not concerned with this issue in this thesis.

We studied Z (t) for various temperatures for a fixed density. These results are
equivalent to those fixing the temperature and changing the density. The rationale
for this assertion is that at a given temperature, a system can change from a liquid
state to a solid state by increasing its density just as one can obtain a solid state
by lowering its temperature at a fixed density. Figure 6.4 is the VAF of 2K = 0.9
system for various temperatures T at a fixed density. 2 (t) at high T (similar to
low density) clearly shows weaker oscillations compared to that at low T (similar to
high density). The weaker oscillatory t-dependence gives a larger positive value for

Z (t) over time which corresponds to a higher diffusion constant.

6.3 Intermediate Scattering Function

Spatial Fourier transform of the density correlation function is more interesting
to study because it can be measured in a laboratory. F(q,t), the spatial Fourier
transform of G(r, t), is called the intermediate scattering function and this quantity
is measured by neutron spin-echo technique [16]. F (q,t) can be separated into

two parts just like G(r,t), namely, the self term (F.(q,t)) and the distinct term
(FJ(Q1t))

F(q,t) = /G(r,t)e“‘"dr
z eaq-(r.(:)-rm(0))

1,171

i
N
= .1- z ”(rm-mo» + i Z eiq-(n(t)—rm(0))

N m N I¢m
= F-(qJ) + F401, t)- (6.22)

123

IIIIIIIIIII'UIIITUIIrUT

1.0

2K = 0.9

Z(t)

 

 

 

LL 1 I l 1 1 l l I l l j l l l l l l l l l l
0 1 2 3 4 5
t (98)

 

Figure 6.4: Velocity auto-correlation functions for various temperatures for 2K =
0.9. From bottom to top curve, temperature is about 250, 400, and 1450 Kelvin.

 

124

The self part describes a tagged particle motion and the distinct part describes the
correlated motion. Figures 6.5-6.7 give the time dependence of F.(q, t) for different
values of q (=0.4, 1.2, and 43“). Each figure has six curves which correspond
to different strengths of the corrugation potential. First, let us understand the q-
dependence of F,(q, t). For a simple homogeneous liquid such as the 2K = 0 system,
F.(q, t) has the following form for small q,

F.”(q,t) = 6'0"". (6.23)

In contrast to that for a homogeneous fluid, F,” (q, t) of a single harmonic oscillator
with frequency can shows an oscillatory behavior and never decays. The amplitude

of oscillation increases as q increases.

_& T 2 -00.
F,”(q,t) = e 353' (1 W". (6.24)
The oscillatory behaviour of F.(q,t) for 2K = 0.9 at short times (Figures 6.5-6.7)
is real and it grows as q increases. This is an effect of the corrugation which results

in a harmonic oscillator-like dynamics of the tagged particle.

Now, let us focus on one of three figures (Figure 6.6) and discuss the time-
dependence of F,(q, t) for different 2K (the three figures 6.5-6.7 show qualitatively
similar time-dependence, so one figure is enough for the purpose of this discussion).
For very short times, these curves show no difference regardless of 2K which can be

understood via a Taylor expansion of F,(q, t) for short times [1].

I
00

8.010 = Z (2 —-)-’2.F.""’(q.0) (6.25)

71:0

,t2 2 t‘ 2
= 1" -w092! +w09wls4!— w09w196!6+ (626)

 

125

 

 

0.6

MM)

0.4

0.2

 

 

 

0.0 ‘

O

 

Figure 6.5: The self part of the intermediate scattering functions at q = 0.424”1 for
different strengths of the corrugation potential. From the top curve, 2K = 0, 0.15,
0.36, 0.55, 0.75, and 0.9.

 

F.(q.t)

1.0

0.8

0.6

0.4

0.2

0.0

126

 

 

I Ifi I l I T I I l I I I I I I I I I I I I I
q
' -1 2A" '
. (1" - .
I .
~O
.; .. .
I- ‘K \ -4
)- }.-. \ -4
—‘\ ‘. \ —
. (\2 ~ .
.s -‘
’ 1‘ . s ‘
"’ 3‘ \. s q
)- 5‘ . \ ‘ -
s...
— ".\\ .\ ....... -
- s. x - ------ -
.P \ \ .°‘
b I.\ \ I
... .‘
- '.\ \ ' nl
\ \
h \ \\ 0‘ ..... II
— \ '.. _
\ 's.
[- \\\ .._~~
\ \
‘ \ x“ "
)- ‘ “~ .1
o:‘ ~~~~~
b ’0..: ..... :
1 1 1 1 I 1 1 1 1 I 1 1 1 1 If..1"1::f°§3l‘w

 

 

2 4 6 8 10
t (p8) -

Figure 6.6: same as Figure 6.5 except q = 1.22214.

 

F.(q.t)

127

 

1.0 IIII'IIIIITIIIIIIIIIIIII

L

q=4A"1

0.8

I'
l

0.6

Ijjr'

0.4

 

lLllllllLlllllLllllj

J I I I I

0.2

   

'2‘?

0‘ \.

I I I I
‘
s
s
I
O
I

 
   

 

 

.ati

0.0

Figure 6.7: same as Figure 6.5 except q = 4A".

128

 

where
w?» = q’%. (6.27)
wi“. = w3.d+03, (6-28)
1
2 _ 2
no - ,kaTam. (6.29)

where d is the dimension of the system and F is the total force acting on the particle.
Eq. (6.27) clearly states that the short time behavior of F.(q,t) is independent of
the interparticle potential. Actually, this fact is easily understood because for a very
short time the tagged particle cannot feel the existence of other particles and the
effect of single particle potential, so it is moving like a free particle described by its
average kinetic energy, or equivalently by the temperature (T) of the system through
equipartition theorem. Compare the coefficient of the t2 term which involves only
IcBT with that of the VAF (Eq. (6.21)) where the coefficient of the t2 term already
involves the mean square acceleration. This implies that velocities of particles reach

equilibrium faster than their coordinates.

At intermediate time scales, interparticle interactions start to contribute to the
correlation function. We can clearly see the effect of the interaction potential on
F,(q,t). Also, as the strength of the corrugation potential increases, the function
becomes oscillatory and also decays slowly. As will be discussed later, the oscillatory

motions give rise to broad peaks at finite frequencies in the dynamic structure factor.

The decay of the F,(q,t) with time slows down as the strength (2K) of the
corrugation potential increases. Recall that the F.(q,t) of a harmonic oscillator
[14] with angular frequency «20 does not decay (see Eq. (6.24)). The system with
corrugated potential can be approximated by a harmonic oscillator model except
that the barrier height for this system is finite unlike the harmonic oscillator. Once

a particle gets trapped near the minima of the corrugation potential, it takes some

 

129

time for the particle to overcome the barrier of the corrugation potential. This life
time depends on how high the barrier is compared with the energy of the particle.
The intermediate scattering function which relaxes rather slowly in the presence
of a strong corrugation potential does not decay completely within our total MD
run time. Therefore, finite MD run time corrections should be considered seriously
in calculations of the Fourier transform of time correlation functions. This will be
discussed in detail in the following section which deals with the incoherent dynamic

structure factor.

6.4 Incoherent Dynamic Structure Factor

In general .the dynamic structure factor 5 (q, V) which is time-Fourier transform of

F(q, t) (see Eq. (6.22)) is related to the differential neutron scattering cross section

 

 

 

 

 

 

[1, 7]-
d’a' J’a dfa
dfldE' ‘ dfldE' .9. + dfldE’ , (6'30)
016’ aJc’

- ZEN-9%") + “kNS-(qw), (631)

where
S'(q,1/)=/oo F(q,t)e"2”‘“dt, (6.32)
S.(q,u) = /°° F,(q,t)e"2""‘dt. (6.33)

Here, subscripts coh and s are abbreviations for the coherent and the incoherent
(or self) respectively, E’ is the final neutron energy, I: and k’ are the initial and the
final neutron wave numbers, N is the total number of particles in the target. By
varying the isotopic composition of the sample or by using polarized neutrons, it is

possible to measure the incoherent and the coherent cross sections separately which

 

130

describe single particle and collective dynamics, respectively. The accessible wave
numbers in inelastic neutron scattering experiments lie typically between 0.1 and 15

A“, which is the same range of q that is usually studied in MD simulations.

Before we discuss our simulation results for various 2K, let us briefly review the
experimental results. Zabel and coworkers [7, 8] have carried out neutron scattering
(both quasi-elastic and inelastic) experiments in stage-2 alkali metal GICs in the
fluid phase of the intercalants for a range of q values between 0.75 and 4A". The
three most significant observations in their experiments are: (1) a diffusive central
peak indicating liquid-like dynamics, (2) a broad finite frequency peak (V ~ 1 THz)
whose intensity grows with the scattering wave number q, and (3) an extremely
narrow resolution limited central peak whose width is considerably smaller than that
of the diffusive central peak. Zabel et al. [7, 8] associated the narrow quasi-elastic
central peak with a ‘solid-like’ dynamics, characteristics of a domain-wall solid above
the melting temperature. Fan et al. [9] tried to understand the origin of this ‘solid-
like’ response by carrying out MD simulations for R0024 using a model consisting
of a repulsive screened Coulomb system in the presence of a graphite corrugation
potential. From the long-time behaviour of the self intermediate scattering function
F.(q,t), Fan et al. argued that the resolution limited narrow central peak was
caused by a ‘solid-like’ structure with a life time ~ 26ps. However, the microscopic
nature of this structure, sensitivity of this life time to the MD run time and to the
strength of the corrugation was not addressed in their study. Thus, what causes this
extremely narrow peak in the dynamic structure factor still remains an intriguing

question.

To understand the physical origins of (2) and (3) above, and how these and the

diffusion rate depend on the corrugation strength and the scattering wave number

 

~-—<'

131

q, we have carried out extensive MD simulations in a similar repulsive screened
Coulomb system but by changing the strength of the corrugation systematically
and taking into account the effects of finite MD run times in our analysis. To our
knowledge, this is the first study of the spectral evolution of the dynamic structure
factor as one goes from a 2D homogeneous fluid (uncorrugated system) to a 2D
latticefluid (strongly corrugated system), the corrugation potential dominating the
dynamic behavior of the latter system. From our q dependent studies, we show that
the observations (2) and (3) are intimately connected and are important signatures
of a lattice-fluid. We also show that the anomalously narrow central peak is a
characteristic of a single particle dynamics in the presence of a corrugation potential.
In the absence of corrugation, as we approach the solid from the liquid phase, some
of the fluid particles get trapped in a cage-potential but this does not give rise
to a narrow central peak. In addition, we analyze and explain how the diffusive
central peak of the incoherent dynamic structure factor narrows as one goes from a
homogeneous fluid to a lattice-fluid and how the intensity of finite frequency peak

changes with q.

From S,(q, V) in the liquid state, one can extract the diffusion constant in the
hydrodynamic limit [1, 3, 14]. However, neutron scattering experiment is not a
direct way to obtain the diffusion constant if the sample is a dominant coherent-
scatterer like Rb in R0024 system. In this case, one uses an approximation [17] which
relates S (q, V) to S.(q, V) through the relation S,(q, V) = S (q, V) where 'q’ = q/m
and S (q) is the static structure factor. This approximation is reasonable only when
vibrational motions (Brillouin peaks) are well separated from the central peak of
S.(q, V). MD simulation is a good method to avoid this approximation because

using MD it is possible to calculate S.(q, V) directly. In this section, we will focus

132

on the S.(q, V) to study the single particle dynamics. Different systems such as a
classical harmonic oscillator, a simple fluid (2K = 0), and systems under the effect
of the corrugation potential (host layers or substrates) are considered and results
will be compared with one another to understand the effects of the corrugation on

the single particle dynamics in 2D.
6.4.1 Classical Harmonic Oscillator

As a guide to understand the liquid dynamic structure factor in the presence of
external corrugation potential and how it evolves from the dynamics of a solid,
let us look at the dynamic response of a classical harmonic oscillator. When the
temperature of a system is low enough, we can use harmonic approximations to study
the dynamic response of the system. In the case of a single harmonic oscillator [14]
with a vibrational frequency 00, one has ’

Sn(q.u)=exp(-y) f 6(u-nuo)1.(y)= i 6(v—nuo)A.<y). (6.34)

732-” n=—oo

where y = chqu/m(21rVo)2 and I,,(y) is the modified Bessel function of the first
kind. The frequency V0 is given by JEFF/211', where C is the second derivative
of the potential at its minimum. If we apply this to the system with corrugation
potential, V0 = 0.982 x J27 (THz). EQLL—‘Jl, Ao(y) = l for y = 0 and decreases
with 3] (see Figure 6.8) whereas fQLflil A1(y) = 0 for y = 0 and increases with y
(see Figure 6.9), namely, for a given 2K (i.e. V0) the spectral weight of the central
peak (finite frequency peak) decreases (increases) as q increases. Furthermore, if
we compare Figures 6.8 and 6.9 for a given y (i.e. T and q), we see that An(y)
decreases as n (i.e. V) increases. For small q, An(y) decreases very quickly with

n and for large q, An(y) decreases with a Gaussian shape as a function of n (or

 

133

 

 

 

 

 

1.0 I 1 I I I 1 I I I I I I I I m I I
; . 1 7
7- o D u
)- O DD+++ -
0.8 '- 60 no ++ "‘
_ o c1:1 ++ :
r 0 DD +I+ '1
- 0 D '0“. '4
0.6 — 9 on *+ d
f; - 0 DD ++ +
v ’ 9 D +* '1
4‘1, : 00 can +++ -
0.4 ’— 090 D'.. +++++
i x 2K = 0.9 3..o ’*+
g 8 2K = 0.75 °°.,° ' -.
0.2 -_— + 2K = 0.55 °°°°o¢°°°°
. 0 2K = 0.36
; 0 2x = 0.15
0.0 l l l l l j l I l _l_l l 1 l l J l l
0 1 2 3 4
‘1

Figure 6.8: Spectral weight of the incoherent dynamic structure factor of a harmonic

oscillator at V = 0 for various strengths of the corrugation potential. y in y-axis is
11:)3Tq2/m(21rVo)2 where T = 250K and m is mass for Rb.

134

 

 

     
 
  

 

 

0.30 bT I I I I I I TI I I I I I I I I 1 r I d
. J
0.25 T 1
- 4
: o°°°°°o
0.20 — .° on
C ° 0
A " o on +
Q0.15 F" 6° DD +++x x
< : 0 OD ++ xxx):
.. D + xxx
0.10 "' o D '3 xx
)- D + xx
; n ++x:x 3‘ 2K = 0.75
. ° ° + xx + 2K = 0 55
_ a + Xx '
0'05 ; 0°00 + "" a 2K = 0.36
-. ’ 02K=0w
0.00 b I I l I I I I l I L l I l I I I
0 1 f 2 3 4

‘1

Figure 6.9: Spectral weight of the incoherent dynamic structure factor of a harmonic
oscillator at V = V0 = 0.982 x \/ 2K for various strengths of the corrugation potential.
y in y-axis is kBqu/m(27rVo)2 where T = 250K and m is mass of Rb.

135

V) although the peaks occur at discrete values of V = nVo, where n is an integer.

Therefore, it is enough to consider only a few low frequency peaks of S g(q, V).

6.4.2 Corrugation-free System

Let us now discuss the dynamics of a corrugation-free (2K = 0) fluid. In Figure 6.10,
we give .S',(q, V) vs. V for different values of q starting from q = 0.1 to q = 4 (in units
of A") at T = 254K. The liquid to solid transition temperature for this system is

about 200K [19]; thus we are not too far above the transition temperature.

In the small q and V limit, the incoherent dynamic structure factor has the
well-known Lorentzian shape appropriate for a liquid. The Lorentzian form for
S.(q, V) is obtained by using the Gaussian approximation for the time dependent
pair distribution function G.(r,t). The spatial Fourier transform of G,(r,t) then
becomes F,(q, t) = e’an‘. Since time correlations in MD run are reliably known up
to t,, where 2t, is the total MD run time, we must consider this finite time range in
calculating the time Fourier transform of G,(r,t). Taking the Fourier transform in

the time interval (0, t,), we obtain

2Dq2
(.02 + (Dq2)2

249

52 + (0.2)28‘D"'" 81.157.76.35)

 

[1 — e'qut' cos Lair] +

 

5.01.”) =

where w = 277V. When t, —+ 00, Eq. (6.35) recovers the well-known Lorentzian

central peak,

2Dq2
S.(q, V) = w” + (Dq2)2' (6.36)

 

Therefore, terms with cosine and sine in Eq. (6.35) are due tofinite t... If we take

V = 0 in Eq. (6.35), then

2 2
_ _ -D r
S.(q,l/ — 0) — [I — 6 q t ] X D—qi' (6.37)

136

 

q
l

IIIIIIIIIIIII'IIIIIIIIIIIII

t"
o

2K=0

.°
a:

S.(q.V)

IIIIIIIIILLIIIIIIIIIIIIIIL

 

 

0.0 * “ “ ‘
0.0 0.5 1.0 1.5 2.0 2.5 3.0
V (THz)

Figure 6.10: Incoherent dynamic structure factor vs. frequency for a system without
corrugation (2K = 0) in the liquid phase at 254K. Various curves are for different
q (in .21"). q increases in steps of 0.1 from 0.1 to 2.0, excepting the broad curve
which is for q = 4 showing Gaussian-like response appropriate for a free particle.

 

137

When Dq’t, is less than or comparable to 1, it is important to take into account the
effects of finite MD run time in obtaining the diffusion constant D from S.(q, V).
For the system 2K = 0, this correction term is negligible. However, the diffu-
sion constant D decreases as 2K increases, therefore the correction factor becomes

significant for large values of 2K.

By fitting S.(q, V) calculated with coordinates obtained from MD simulation to
Eq. (6.35), we can obtain diffusion constants. We repeat this fitting for various
values of q from 0.1 to 1A”. Figure 6.11 gives the q-dependence of D(q) for various
strengths of the corrugation potential. The diffusion constant D is determined by
extrapolating D(q) towards q —» 0. The value of D agrees well with that obtained
from the relation ((Ar)’) = 4Dt for large t, if we do not consider the value of
D(q = 0.1) seriously. Because q = 0.1 corresponds to a length scale which is about
half of the MD cell size (recall that periodic boundary condition is used), finite size

effect is more serious when D is large. ‘

Continuing with the discussions for a corrugation-free fluid, for intermediate
values of q, we see a weak broad peak growing at finite frequency near 0.9THz
(see Figure 6.10). The origin of this finite frequency peak can be ascribed to the
dynamics of a particle trapped in the field of its nearest neighbors, the so called
cage dynamics. As far as we know this is the first clear MD observation of the finite

frequency cage dynamics in a 2D homogeneous fluid.

For large values of q, the central peak becomes a Gaussian which characterizes
a “free” particle dynamics. Therefore, there is no qualitative distinction between
diffusional and vibrational motion, both occurring on the same time scale. This
dynamic structure factor at q = 4A" is Gaussian-like as expected for a homogeneous

fluid at large values of q and there is no evidence of a narrow central peak as seen

 

D(q)

Figure 6.11: Diffusion coefficients from S.(q,V) as functions of q. D(q —v 0) ap-
proaches the value of D obtained from other methods (see Table 6.1). The symbols
x and D at q = 0 are the values of D obtained from mean square displacements

0.4

0.2

0.0

0.0

138

 

 

IITIlrIIIIIIIIlIIIITIIII

.1

 

 

IIIIlIIIIlIIIIlIJILiIIII

0.2 0.4 0.5 0.8 1.0
q (A")

and the Fourier transform of the VAF, respectively.

 

139

in Rng; neutron scattering experiments [20].

6.4.3 Effect of a Strong Corrugation Potential

Now let us discuss the effect of a strong corrugation potential (2K = 0.9) on the
dynamics of the fluid. Figure 6.12 gives S,(q, V) for T=248K. To get a physical feel-
ing for the strength of this corrugation, we find that the barrier height A associated
with only the single particle potential is 1080K. For this system the transition tem-
perature from the liquid to the solid phase is about 220K [19]. Thus we are probing
the dynamics of a highly correlated lattice-fluid since the intercalants spend con-
siderable amount of time near the bottom of the single potential well. As we see
in Figure 6.12, the spectra for 2K = 0.9 is drastically different from those for the
2K = 0 case (see Figure 6.10). The narrow central peak is indicative of Rb dif-
fusion occurring on a time scale much slower than that of the lattice vibrations
corresponding to the finite frequency peak. We find that as q increases, the height
of the narrow central peak decreases while it broadens, and in addition the strength
of the finite frequency peak increases. The diffusive Lorentzian central peak for
small q narrows considerably in the presence of corrugation indicating a dramatic
reduction in D. Fitting the spectrum to a Lorentzian with finite MD run time
correction (Eq. (6.35)), which plays a very important role in this case, we obtain
D = 0.03 x 10"cm’/sec. This value of D agrees well with a calculation of D from
the relation ((Ar)’) = 4Dt for large t.

When we increase q, a broad finite frequency peak (actually a doublet) grows
near V ~ 1THz. For large q 2 4 i.e. when qa >> 1, unlike the homogeneous
fluid phase, the spectral function for the corrugated fluid deviates drastically from

a Gaussian form. It has a narrow central peak sitting on the top of a broad central

140

 

 

 

 

S.(q.V)

 

 

Figure 6.12: Incoherent dynamic structure factor in the liquid state (at 248K) vs.
frequency for 2K = 0.9. Various curves are for different q’s (in A") which are the
same as in Figure 6.10. The uppermost curve (q = 4) is no longer a simple Gaussian.
The inset gives intensities of the central peak at q = 4 as a function of the MD run
time; 2t, is the total MD run time.

141

peak and in addition there is a broad finite frequency peak. These observations are
in agreement with the inelastic neutron scattering measurements in RbCu by Zabel

et al. [8, 7, 21] and earlier MD simulation studies by Fan et al. [9] (see the following

paragraph for more discussion).

For 0.75 < q < 2.2, Zabel et al. [7] could not fit their central peak to a single
Lorentzian and suggested that the central peak actually consisted of two Lorentzians,
one representing the usual liquid diffusion and the other (extremely narrow and
resolution limited) corresponding to a ‘solid—like’ peak. For q = 4, the narrow solid-
like central peak was still observed on the top of a broad liquid-like central peak [8].
Fan et al. [9] in their MD simulation study at 300K, find that S,(q, V) has a very
narrow central peak whose life time T is about 26 ps which is about half of their
maximum MD run time (50ps). As discussed before, correlation functions in MD
simulations can be reliably calculated up to about half the maximum MD run time.
Therefore one has to take the value of 1' = 26ps with some caution. Furthermore,
the precise physical origin of this narrow peak, whether due to collective or single-

particle effects, was not probed by Fan et al..

To answer the questions posed above and to determine the sensitivity of 1' to the
MD run time, we have analyzed S,(q, V) in considerable detail and also by increas-
ing the maximum MD run time up to ~ 140ps. We also find that for large q [21],
S.(q, V) shows a very narrow central peak suggesting that there is an additional con-
tribution to the dynamic central peak beyond the usual liquid dynamics. Assuming
the existence of a ‘solid-like’ component with life time 1' and allowing for the finite
MD run time (2t,) correction, we fit S,(q = 4, V = 0) to Cl[1 — e"'/'] + 02, where
the contribution 02 comes from the usual fluid dynamics. Finite time correction to

the second term is negligible because correlation function for a fluid dies out rather

142

rapidly for these large values of q. In Figure 6.12 (inset), we plot S.(4,0) vs. t,
(=18.53, 25.31, 33.8, 40.59, 46.53, 55.86, 60.95, 67.74, and 70.57 ps). The fact that
the slope of this curve decreases with t, suggests that the narrow central peak is not
a 6—function. For a 6-function time response with finite MD run time correction,
the slope would have been linear. We find r=33.4 ps which is slightly larger than
26 ps found by Fan et al. [9]. One reason for this difference is perhaps the lower
temperature in our simulation than theirs and the other is the limitation in their

(Fan et al.) MD run time.

6.4.4 Evolution of S,(q, V) with Strength of Corrugation

To see how the spectral function .S'.(q, V) evolves from 2K = 0 (homogeneous fluid) to
2K = 0.9 (lattice-fluid), we plot in Figure 6.13 S.(q, V) for q = 4021’1 and 0.4.3"1
for different values of 2K near 250K. For q = 4A“, i.e. in a short length scale
(qa >> 1), the complex spectrum corresponding to 2K = 0.9 evolves to a Gaussian
for 2K = 0 by modifying all the spectral characteristics. The finite frequency peak
softens and broadens and the central peak broadens. The central and finite frequency
peaks are no longer distinguishable from each other. For 2K = 0.15 (A = 180K)
corresponding to the condition kBT > A (see Table 6.1 in page 117), the spectrum
is close to a Gaussian but for 2K = 0.36 (A = 432K) and larger i.e. kBT < A,
a narrow central peak appears and its intensity grows with increasing strength of
the corrugation and simultaneously a broad finite frequency peak starts to grow.
When A is increased still further (e.g. up to the lattice gas limit), we expect the
central and finite frequency peaks to be narrow and approach a 6-function. As we
increase the temperature of a system with a nonzero 2K until the thermal energy of
a tagged particle becomes larger than the barrier height A, we find that the intensity

of the finite frequency peak decreases and finally becomes indistinguishable with the

143

 

 

 

 

S.(q=4,l/)

 

 

 

 

' 0.0 0.5 1.0 1.5 2.0
V (THz)

Figure 6.13: Incoherent dynamic structure factor vs. frequency for various strengths
of the corrugation near 250K for q = 4A". This finite frequency dynamics is due
to the trapping of particles near the minima of the corrugation potential. At low
frequency, we can clearly see the evolution of a sharp peak when the barrier height
A is higher than IcBT. Inset: Incoherent dynamic structure factor vs. frequency
for various strengths of the corrugation near 250K at q = 0.4231“. This shows how
diffusive motion slows down in the presence of the corrugation potential. We see a
dramatic increase in the central peak intensity as 2K increases from 0 to 0.9.

144

broadening central peak and the overall shape approaches a Gaussian.

To probe the nature and origin of the resolution-limited central peak more pre-
cisely, we can go to a small wave vector. This is a distinct advantage of MD simula-
tion over experiment where the coherent Rb scatterers preclude the probing of self
correlation functions in long length scales. For q = 0.4 (inset in Fig. 6.13), the entire
spectrum consists of the central peak with vanishing weight in the finite-frequency
peak. It can be fitted to a single Lorentzian with finite MD run time correction;
its width narrows as we increase 2K . This corresponds to the slowing down of the
diffusion in the presence of the corrugation potential and can be thought of as an
effective mass enhancement of the diffusing intercalants caused by the corrugation
potential. Within the accuracy of our calculation, there is no evidence of the very
narrow ‘solid-like’ component due to corrugation potential in this long length scale.
Thus a proper physical understanding of the origin of the ‘solid-like’ component

must address this strong q-dependence.

A simple way to understand the origin of the narrow quasi-static central peak
and the q—dependent intensity of both this and the finite frequency peak is to recall
the dynamic response of a classical harmonic oscillator discussed in Sec. 6.4.1. The
spectral weight of the finite frequency peak SH(q, V = V0) for a harmonic oscillator
increases as q increases. Therefore, we understand why the finite frequency peak of
2K = 0.9 system increases as q increases. On the other hand, the spectral weight of
the central peak, SH(q, 0), decreases as q increases. For large q, even if the spectral
weight of the central peak A0 decreases with q, we clearly see the narrow cental
peak on top of a broad central peak. Here, we must comparethe spectral weight
of the narrow component with that of the fluid component of the central peak. We

argue that there exists a crossover length scale characterized by a wave number

145

q‘ = J17??? where 1' is the average trapping time and D is the liquid diffusion
rate, such that for q > q‘, the narrow solid-like component dominates the central
peak whereas for q < q‘, the diffusive Lorentzian component is dominant. For the
2K = 0.9 system, we find q‘ z 1 which explains why for q = 0.4 we only see a single
liquid-like difiusive peak whereas for q = 4 we see that the solid-like component
dominates the spectrum. Precisely speaking, S.(q = 0.4, V) is fitted with Eq. (6.35)
within about 5%. If we treat this 5% difference as due to the solid like response, the
Lorentzian central peak is about 20 times bigger than the solid-like peak. On the
other hand, for q = 4, S.(q = 4, V = 0) = C; [1 - e“'/"] + Cg, where the contribution
Cg comes from the usual fluid dynamics. We find that at 250K, C1 = 10.2 and
C; = 0.84. Intensity ratio of the solid-like component and the liquid component is
CEO - e""/’)/Cg = 10.7. Therefore, in the incoherent dynamic structure factor at
q = 4 for the system with 2K = 0.9 at 250K, the solid-like component is about 10

times stronger than the liquid-like component.

6.5 Coherent Correlation Functions

A liquid is a dense medium where fluctuations occur continuously and spontaneously.
Local perturbations from these fluctuations disturb the equilibrium state of the sys-
tem. During its relaxation, each degree of freedom (or mode) of the system returns
to its equilibrium. Those modes whose characteristic decay time is long compared
with the molecular interaction time are the collective modes. These processes involve
many particles and their relaxation time is proportional to the square of their char-
acteristic wavelength which is large compared to the intermolecular distance (the life
time of fluctuation is determined by a thermal diffusivity). The description of these

long-lived collective phenomena constitutes the object of linear hydrodynamics.

ii
,,.

146

To see the effect of corrugation on the collective dynamics in a 2D fluid, we con-
sider the coherent intermediate scattering function F(q, t) and its Fourier transform
which is called the coherent dynamic structure factor S(q, V). F(q, t) and S(q, V)
give a detail description of the collective modes in a fluid. First let us look at the

short time behavior of the coherent intermediate scattering function,

00 t2n
= (2n) .
F ((1.1) E0 (2n)!F ((1.0) (6 33)
_ F 2‘2 2 2t4 2 4‘6 6 39
- (mm-«2.5,+wow.;,-,-wow.-6—,+-~. (- )

where F(q, 0) = S(q) is the static structure factor and it gives us information about
the interparticle potentials indirectly. Here, 093 is equal to 103, in Eq. (6.25) which
can be shown either by using the equation of continuity (mass conservation) or by
simply calculating Fm(q,0) directly. Figures 6.14-6.16 give the 2K dependence of
the coherent intermediate scattering function for q = 0.4,1.2, and 4191'1 near 250K.
The well defined oscillatory motion for q = 0.4/“1’l (Figure 6.14) is a signature of the
collective mode (Brillouin peaks) which appears in the coherent dynamic structure
factor S (q, V). As we increase the strength of the corrugation 2K from 0 to 0.9,
F(q,t) shows shorter period which leads to a shift of Brillouin peak position to
higher frequency. For q = 1.2.?1'1 (Figure 6.15) which is very close to the value of
q corresponding to the maximum of S (q), the large amplitude and the slow decay
of F (q,t) is related to the strong spatial correlations that exist at this wavelength.
However, the oscillatory behavior is nearly absent excepting for very large values
of 2K . Figure 6.16 is quite similar to Figure 6.7 (which describes the incoherent
dynamics) indicating that for q ~ 4271-1, namely for short wavelengths, the coherent
and the incoherent density correlation functions are nearly identical telling that

the self term in Eq. (6.22) dominates the pair terms in the intermediate scattering

147

 

0.02IIIT'IIIIIIIII‘IWIII'TjII
1

q=0.4A'1 '

0.01 -- 5'; -

 

F(q.t)

  

 

0.00

r

J I
.I.”U""-z .e......
000‘.... .
so...sooooo..l.o.o.?....
-
a
"‘\

00".
s

I

r
I

-0.01 43.3 .—

IIIIlIIIIlIIIIlIIkILJIIL

0 2 - 4 6 8 10
t (PS)

 

 

 

Figure 6.14: The intermediate scattering function at q = 0.4/'4'1 near 250K. From
top (solid) curve, 2K=0.9, 0.55, and 0 (dotted curve).

MM)

148

 

   
 

4 I I I I I 1 I I T ‘l I I I I I I I I 1 I I IT I
- t
b. 1
II: q= 102A
\.
3 - a
Q
" \
f.
F O 0‘.
p '. '. '-
. 0‘
P ' u
2 - _
7- ". -I
P ‘. . II
\
b . .‘. I
1- ."o.. q
.....
1 L- ----- ...—.
[- d
D I
#- ... u
0 _ 00.......:
I I I I I I I I I I IA LI. 1 I I I I l I I L I ‘

 

 

 

0 Z 4 6 8 10
t (PS)

Figure 6.15: Same as Figure 6.14 except q = 1,2}1-1

F(q.t)

 

149

 

IIIIIIIIIIIIjIlIrWIITIII

q=4A'1

I

IIIIIILIIIIJIIIIIIIIIII

IILJ

 

 

O 0 1.21.1 LI 1 1 1 1- 1.4 1 1.1.1.?1'1: °'1"l'"f'J"
O

0 2 4 6 6 10
t (PS)

 

Figure 6.16: Same as Figure 6.14 except q = 421-1

150

function.

Now, let us look at the coherent dynamic structure factor 5 (q, V) for different
strengths of the corrugation potential (see Figures 6.17-6.19). At q = 0.4121"1 , we see
the Rayleigh (central) peak and the Brillouin (finite frequency) peak. As we increase
the strength of the corrugation, the Brillouin peak moves to a higher frequency and
the peak intensity decreases. Simultaneously the central peak narrows and grows in
intensity. The increasing frequency of the Brillouin peak indicates that the sound
velocity increases in the presence of corrugation. At q = 1.2131“, the strength of the
Brillouin peak is negligible for all values of 2K, but the narrowing feature of the
central peak persists. As we have mentioned before, at q = 421-1, the incoherent and
the coherent intermediate scattering functions are almost identical. Correspondingly
S (q = 4, V) is practically same as S,(q = 4, V) given in Figure 6.13. The reason that
S (q, V) is more noisy than S,(q, V) is simply because the angular averages were done
exactly for F.(q,t) (analytic integration) whereas for F(q, t), they were carried out

numerical discrete summation.

Finally, let us study the sound velocity of the system by looking at the q depen-
dence of the Brillouin peak which is the finite frequency peak in Figure 6.17. We
can obtain more accurate information about the Brillouin peaks through the longi-
tudinal current correlation function J1(q, V) which is directly connected to S (q, V).

The relation is as follows:

F(q,t) gamma» (6.40)

7.84.0 = fiwmwa». (6.41)

where

nq(t) = Z 69'le (6.42)
I

151

 

 

0.08 h I I I I If I rfi I I I T I I I I 4

. _ _1 .

[ q—0.4A '1

0.06 — '-
3 .
3 0 04 -
”J I

  
     

0.02

l

l I I I

 

 

 

Figure 6.17: Dynamic structure factor for 2K = 0, 0.55, and 0.9 from left to right
peak at finite frequency at q = 0.424”.

152

 

IIFIIIIIIIIIIIIIIIIIIIII]ITII-I

125
q= 1.2A'1

100 _—

q
01
[III

 

0|
0
T'III

 

IIIIIIIIIIIIIIIJIIIIIIIIIL

 

 

:\
25—
’..
I
b
o...
h d

 

 

'111111111'111111111I_1111I1111

0.00 0.05 0.10 0.15 0.20 0.25 0.30
V (THz)

 

Figure 6.18: Same as Figure 6.17 except q = 1,2}1-1

153

 

InorfiIIIlrwvuluvrtlru

i
: q=4.0.4'1

-

0.8

)-
[-

I
1
I
I
4
1
I
4
1
1

0.6 i-

 

S((LV)

IIIJJIIIIJIIIIIIIIIIII

 

 

f

 

0.01111I1111l1111111
0.0 0.5 1.0 1.5

' V (THz)

 

Figure 6.19: Same as Figure 6.17 except q =.- 421-1

154
13190) = sza(t)€’q'n'(t’o (6.43)
1

By combining the two equations

62

5721“ (7.1) = -q’Jz(q.t) (6-44)
gum) = iq . 14(4). (645)

we obtain
Jz(q. V) = (211023 («1. V)/ ‘12- (6-46)

This relation is exact and follows from the number conservation through the equa-
tion of continuity (Eq (6.45)). While this relation implies that the two correlation
functions have the same information, certain spectral properties are more easily dis-
cerned in one function than in the other. In fact, it is more appropriate to discuss
the collective modes in a liquid in terms of Jg(q, V) rather than .S' (q, V). The factor
V2 in Eq. (6.46) effectively suppresses that the central component (Rayleigh compo-
nent) in S (q, V) which does not show up in J;(q, V). If the Brillouin peak in S (q, V) is
sharp, it will appear also as a sharp peak in J;(q, V) at the same frequency. However,
if the peak in S (q, V) is broad, it can appear considerably distorted in J1(q, V) and

the peak position can be shifted to a higher frequency.

Figures 6.20 and 6.21 give the longitudinal current correlation functions as a
function of q (< 0.6.21") for 2K =0 and 0.9. We can see how the peak position
moves to a higher frequency as q increases. For the 2K = 0.9, some shoulders
come from the large central peak and should not be included in the collective mode
response. The dispersion of the peak frequency vs. q gives us information about
the sound velocity of the system in the hydrodynamic regime. [In other words, the
peak position in J1(q, V) can be regarded as the frequency of the sound waves at that

wave number q. Therefore, we plot in Figure 6.22 the q-dependence of the frequency

155 l

 

0.25IIIII’IIIIIIIIIIIIIII‘I'IIII [

I

[ 2K=O
. 0.20

 

0.15

 

Jl(q’V)

0.10

 

 

 

0.05

 

 

 

 

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

 

TIIIIIITITIIIIIIIIIIIII

   

0.00
2.5

Figure 6.20: Longitudinal current correlation functions are shown as a function of
q (from left q = 0.1, 0.2, 0.3, 0.4, 0.5, and 0.6 in A“) for system 2K = 0.

156

 

0.12 IIIIIIII'IIIIIIITIfiIII

2K=0.9

    

0.10

0.08

0.06

J1(<1J’)

 

0.04

0.02

 

IIIIlIIIIlIIIIlIIlIIIIJIIIILI

 

 

2.5

Figure 6.21: Longitudinal current correlation functions are shown as a function of
q (from left q = 0.1, 0.2, 0.3, 0.4, 0.5, and 0.6 in A") for the system 2K = 0.9

 

157

at maximum intensity of J;(q,t). We obtain a sound velocity 1.7 x 105cm/sec for
the system 2K = 0 from the linear dependence of the frequency on q. However, for
large 2K, we find that the collective mode frequency appears to approach a finite
value as q goes to 0. For a liquid the collective mode frequency must be zero for
q ~ 0. Therefore for large 2K there is a rapid increase in V, as q increase from 0,
thus giving rise to a high sound velocity. This is because most of the particles are
trapped inside potential well and it costs a large energy to excite long wavelength
collective modes. But from the finite size simulation, we are unable to get a more
quantitative description of the sound wave dispersion in the small q limit in the
presence of a strong corrugation potential. When we increase the temperature such
that the particles can diffuse easily, we obtain the usual low sound wave frequency.
For example, for the 2K = 0.9 system at high temperatures (1450K) which is above
the potential barriers, the effect of the corrugation potential can be neglected and

the frequency goes to 0 linearly as q —+ 0 as in the corrugation free system (see

Figure 6.22).

6.6 Summary and Discussion

We have studied the effects of the corrugation potential on the dynamic properties of
a two-dimensional fluid. The model system we have studied can represent not only
ionic overlayers on corrugated substrates but also intercalated systems by choosing
appropriate parameters. Neutron scattering studies of the dynamic structure factor
at large wave numbers from a two-dimensional (2D) fluid on a corrugated substrate
show an extremely narrow resolution limited central peak on the top of a fluid central
peak, and a broad finite frequency peak. Using molecular dynamics simulations we

reproduce all these features which are understood as general characteristics of the

*-

158

 

1.5OIIIIIII1IIIIIIIIIFII IT—II

I

1.25

1.00

§ 0.75

0.50

0.25

TIIIITIIIIIIIIIIITIIIIIIIII

 

 

IIIIIIIIIIILIIIIIIIJIIILJIJII

 

IIIIIIIIIIIIIIIIIIIIIII

00
0.0 0.2 _0.4 0.6 0.8 1.0
‘1

Figure 6.22: q-dependence of the frequency at maximum intensity of the Brillouin
peak near 250K except that only the solid circles are at 1450K for the 2K = 0.9
system. From top to bottom 2K = 0.9, 0.75, 0.55, 0.36, 0.15, 0.

 

159

dynamics of a 2D lattice-fluid. Evolution of these spectral features as a function of
the corrugation strength and scattering wave number clearly brings out the physical
origin of the narrow ‘solid-like’ central peak. We have shown how the dynamics of
a homogeneous fluid changes to that of a lattice-fluid where the fluid particles are
trapped for sufficiently long time before undergoing an activated diffusion process.
The corrugation potential which traps the particles with an average trapping time
1' gives rise to the ‘solid-like’ narrow central peak which can be seen clearly at large
wave numbers. In the absence of the corrugation potential i.e. for a homogeneous
fluid, even when we approach the transition to the solid phase we find no evidence

of such a narrow peak at large q.

In addition to the single particle dynamics, we also studied the collective modes in
the liquid state through the coherent dynamic correlation functions. For small q, we
saw both Rayleigh and Brillouin peaks. The Brillouin peak is well defined enough to
study the sound wave dispersion curve. For the corrugation free system, we obtain
its sound velocity from the dispersion curve. As we increase the strength of the
corrugation, the sound velocity of the system increases. We hope that experiments
can be done in the corrugation free system to confirm our results. Also of interest
will be measurements of the sound velocity propagating in the layer direction in

RbC'24 system.

 

Bibliography

[1] J.P. Hansen and LR. McDonald, Theory of Simple Liquids, (Academic Press,
London, 2nd edition, 1986).

[2] M.P. Allen and DJ. Tildesley, Computer Simulation of Liquids, (Oxford Uni-
versity Press, Oxford, 1989).

[3] Jean Pierre Boon and Sidney Yip, Molecular Hydrodynamics, (McGraw Hill Inc.,
1980).

[4] J. Koplik, J .R. Banavar, and J .F. Willemsen, Phys. Rev. Lett. 60, 1282 (1988);
M. Sun and C. Ebner, Phys. Rev. Lett. 69, 3491 (1992).

[5] RA. Thompson and Mark 0. Robbins, Phys. Rev. Lett. 63, 766 (1989) and

references therein.

[6] Graphite Intercalation Compounds, Vol 1: Structure and Dynamics, H. Zabel

and S. A. Solin, eds. Springer Series on Topics in Current Physics, Springer (1990).

[7] H. Zabel, A. Magerl, J. J. Rush, and M. E. Misenheimer, Phys. Rev. B 40, 7616
(1989).

[8] H. Zabel, A. Magerl, A. J. Dianoux, and J. J. Rush, Phys. Rev. Lett. 50, 2094
(1983); H. Zabel, et al., Phys. Rev. Lett. 57, 2041 (1986).

[9] J. D. Fan, George Reiter, and S. C. Moss, Phys. Rev. Lett. 64, 188 (1990).
160

11..

 

161

[10] W.A. Curtin, Phys. Rev. Lett. 59, 1228 (1987); D. Zhu and J.G. Dash, Phys.
Rev. Lett. 60,433 (1988).

[11] MB. Reiter, D.D. Awschalom, and N.W. Schafer, Phys. Rev. Lett. 61, 966,
(1988); NJ. Wilkinson et al., Phys. Rev. Lett. 69, 3535 (1992).

[12] SW. Lovesey Dynamics of Solids and Liquids by Neutron Scattering p.1,
Springer-Verlag, Heidelberg, (1977).

[13] van L. Hove, Phys. Rev. 95, 249, (1954).

[14] W. Marshall and SW. Lovesey, Theory of Thermal Neutron Scattering, (Oxford
University Press, Oxford, 1971).

[15] B.J. Alder and T.E. Wainwright, Phys. Rev. Al, 18 (1970).

[16] M. Davidovic and A.K. Soper, Static and Dynamic Properties of Liquids,
Springer Proceedings 40, (1989).

[17] K. Skéld et al, Phys. Rev. A 6, 1107, (1972).

[18] Hyangsuk Seong, Surajit Sen, Tahir Cagin, and S. D. Mahanti, Phys. Rev.
Rapid Communications B45, 8841 (1992).

[19] Hyangsuk Seong, S. D. Mahanti, Surajit Sen, and Tahir Cagin, Phys. Rev. B
46, 8748 (1992).

[20] W.A. Kamitakahara and H. Zabel, Phys. Rev. B 32, 7817 (1985).

[21] Although Zabel et al. measured coherent scattering structure factor, for q =
4.4“, the coherent and the incoherent structure factors are nearly the same.
Furthermore we find that with increasing temperature the doublet peak structure

becomes a single one as seen in Zabel et al. experiment [7].

Chapter 7

Multiphonon Response in 2D
Corrugated Systems

7 .1 Introduction

As discussed in the introduction and in more detail in chapters 4 and 5, atomic and
molecular overlayers on corrugated substrates and intercalants inside host layers ex-
hibit a rich variety of structures depending on the size of the atoms (or intercalants),
planar density, and the strength and period of the substrate (or host) corrugation
potential [1, 2]. Melting and dynamic properties of some of these physical sys-
tems which exhibit incommensurate ground state structure (or a discommensurate
domain-de structure) have been of considerable experimental and theoretical in-
terest [3, 4]. For melting of an incommensurate solid, the main interest was to
understand the physics underlying the domain and domain-wall melting [3, 5] which
we covered in chapter 5. As regards the dynamics of the incommensurate solid, some
of the fundamental questions related to, for example, the existence of low-frequency
modes and how the in-plane vibrational modes reflect and probe the discommen-
suration structure and the atom-substrate interaction [6, 7, 8]. In this chapter we

address these dynamics questions with the intercalants inside graphite host layers

162

 

163

as our model systems.

Kamitakahara and Zabel (K2) [7, 8] investigated the phonon density of states
(PDOS) associated with the in-plane vibration of the intercalants in a large class
of alkali metal GICs using inelastic neutron scattering experiments with scattering
wave number (3.5 < q < 4.5.4“) for which the incoherent approximation [9] was
found to be justified. Here we focus on the stage-n (n _>_ 2) GICs (in particular
RbCu) where 2-dimensional models have been quite successful [10]. The PDOS for
Rng4 extracted by KZ (to be denoted as g(V), see Figure 7.1) shows two peaks, one
near 1THz and the other near 2 ~ 2.6THz. Even above the melting temperature
(experimentally, Tm ~ 162K) when the intercalants form a 2D lattice-fluid, some
anharmonicity becomes noticeable, but nothing dramatic happens to g(V). The
frequency of the first peak near 1THz was interpreted by K2 to be proportional to
the strength of the Intercalant-Carbon (I-C) in-plane interaction and was therefore
expected to depend critically on the location of intercalants with respect to the
graphite substrate. As pointed out by KZ, this low frequency peak was surprisingly
well defined at low temperature in spite of the fact that many of the intercalants
in the discommensuration domain structure [10, 11] did not sit precisely at the

commensurate sites and therefore’experienced a distribution of LC interactions.

KZ associated the two main peaks with the in-plane transverse (TA) and longitu-
dinal (LA) acoustic modes respectively. They observed that the low-frequency peak
(TA) was still dominant above the melting temperature whereas the high-frequency
peak (LA) broadened remarkably. This is very surprising because even if we consider
the effect of corrugation in the liquid state, we expect the LA modes to be still bet-
ter defined than the TA modes. Furthermore, KZ could fit g(V) obtained from the

lattice dynamics calculations with their experiments in the stage—1 K and Rb GICs

 

164

 

 

 

4 WIIIIIIIIITIII'IIIIfileI
3 000 Liquid
0 000
O
2 o 062,0
00 O
1 0 000m
0 QDQJO 0o
A .0 O
3.0 =
:‘05
4 0 Sand
0
3 O
0
2 o 0 0
0 6136-130
1 o 696’ c?
0 IqQJI1111IJI11oJm111I11
0 1 2 3 4
v(THz)

Figure 7.1: Intercalant in-plane mode phonon density of states from Kamitakahara
and Zabel [7] for RbC24. Top (bottom) panel is in liquid (solid) state.

 

165

[7, 8, 12] only if they used an unscreened Coulomb interaction. They assigned the
high frequency peak to collective mode (plasma oscillation) of unscreened K ‘1’ and
12b" ions. If one assumes that there is no screening in stage-1 GICs, this assump-
tion becomes even more reasonable for the higher stage GICs because of their lower
electronic density. Therefore one should be able to fit even better the 2 ~ 2.6THz
peak in RbCu using the unscreened Coulomb model. However, we calculate this
collective mode frequency to be ~1.35THz thus making the unscreened Coulomb
model inadequate. Also it is difficult to accept the no screening picture in alkali

metal GICs because of the metallic nature of the system.

In their MD simulation studies for Ran using a screened Coulomb model, Fan
et al [13] found two peaks in the self dynamic structure factor S.(q, V). The first
peak was pronounced and centered around 1THz while the second peak was very
weak and centered around 2THz. Fan et al. following KZ’s suggestions, tentatively
assigned the first peak to TA and the second peak to LA mode without studying
how these peaks changed with temperature (by going to the solid state) to confirm

this identification.

In this chapter we propose a completely different explanation of the observed
multipeak structure seen in neutron scattering experiments. We show that although
in general it is difficult to observe multi-phonon response, particularly in the liquid
phase, the substrate corrugation significantly enhances the strength of multiphonon
response. we present the results of our calculation of PDOS of the stage-2 Rb
GIC obtained via both Fourier transform of the velocity auto-correlation function
(VAF) which we will refer to this as gH(V), and also from the incoherent dynamic
structure factor following a procedure that KZ used (defined g(V)). We have used

MD simulations of a 2-dimensional repulsive screened Coulomb system on corrugated

 

166

layers to obtain both gH(V) and g(V). This model, as discussed in earlier chapter
of this thesis, has been very successful in describing the structure [10], melting
[14], and dynamics in the fluid phase [15] of the stage-2 Rb GIC. Our present MD
simulation results for g(V) not only agree with that obtained from the inelastic
neutron scattering experiments [7] but also elucidate the origin of the two-peak
structure. Details about the potentials and MD simulations have been already

discussed in chapter 4.

In order to see how the phonon dynamics depends on the strength of the cor-
rugation and the intercalant density, we will consider three systems; system I with
planar density 0.0318 per A2 in the absence of substrate potentials i.e. no graphite
layers, system 11 with the same planar density as the system I but in the presence
of host layers (or substrates), and system 111 with 0.0311 / A” whose ground state
structure is the periodic domain wall structure discussed that in chapter 4. This
classification is the same as that in chapter 5. The systems II and III are incommen-
surate and commensurate solids, respectively. The strength of corrugation is chosen
to be 2K = 0.9, appropriate for Rng... Comparison of the dynamics of systems with
slightly different densities tells us about the role of atomic defects in the vibrational
excitations of commensurate and incommensurate solids. Layout for this chapter
is following. First, we will discuss the PDOS through lattice dynamics calculations
for systems I and 111. Next the PDOS for the system 11 will be obtained from the
Fourier transform of the VAF and through the incoherent dynamic structure factor.
Finally, the results of a harmonic oscillator model will be presented to illustrate the
wave number (q) dependence of a quantity related to the PDOS but obtained from

the incoherent dynamic structure factor S,(q, V).

167

7 .2 Lattice Dynamics

In this section, we briefly review the theory of the vibrations of a crystal lattice [16].
The atomic displacements are assumed to be small so that the forces are regarded
as linear functions of the displacements. The force constants between atoms are
known as the second derivatives of the interatomic potentials with respect to the
atomic displacements. This is the usual harmonic approximation where the lattice
vibrations are treated as a collection of coupled simple harmonic oscillators. The
Hamiltonian of the system is

H = '3: Mali,“ '1'; E ‘I’aaanseueeeuajm (7-1)

aria: aimfiju

where

Toma,“ = [8’4/6uagxaugjn] (7.2)

0,

<I> being potential energy of the system. For system I, <1) is just the repulsive screened
Coulomb potential whereas for systems 11 and III, the corrugated potential is also
included (see chapter 4) in <1>. um-"s are the displacements of atoms from their
equilibrium positions and the subscripts stand for the ath rectangular component
of the nth atom in the ith unit cell. The equations of motion for the atoms are then
easily found to be

Mnfim'n = - 2 ‘I’as,ie,juusju- (7-3)

3.714

Choosing periodic solutions to Eq. (7.3) (using Bloch’s theorem)

1 -iwt+iqu.- , (7.4)

"air: = _‘/M=nuan(Q)e

we obtain

 

168

wzuanUI) = Z Dafi.nu(Q)uBu(Q)s (7-5)
B»
where
1 - .
Da K = Qa in ' c-IQ’(Ri-R’). 7.6
{3. 14(1) MgMu Reg,- 16. an ( )

B is referred to as the ‘Dynamic Matrix’ and q is the wave vector of the phonon.
For a crystal with r atoms in a unit cell and each atom having d degrees of freedom,
the number of equations to be solved is rd [Eq. (7.5)].

We have obtained the dispersion relation 0) = w,-(q) (j = l, 2, . - - , 2r) by solving
Eq. (7.5) for the systems I and 111 where we have used d = 2 by considering only the
motion of atoms parallel to the plane. Since system 11 has no periodicity (incom-
mensurate solid), it is not possible to use a simple unit cell with a small number of
basis atoms for lattice dynamics calculations. We can form a unit cell with 7-atom
basis for the system 111. Therefore, to directly compare the normal mode frequencies
and PDOS in systems I with those of the system 111, we also choose a 7-atom basis
conventional cell for the system I (see Figure 7.2). Figure 7.3 shows the dispersion
curves for both the systems along the (10) direction of its unit cell. We have 14
branches which come from 7 (=r, number of basis atoms) x 2 (=d, degrees of free-
dom). The figure also shows the linear dependence of the two acoustic modes (LA
and TA) for small q (long wavelength limit) for the system 1. However, the phonon
spectrum for the system III whose ground state is a commensurate structure shows
a gap due to pinning by the corrugation potential [8]. Furthermore the phonon
dispersion becomes very narrow in the presence of the corrugation potential. Once

we know the dispersion curves, we can obtain the PDOS, gH(V), through

9H(V) = 5],— :60 - 16(9)). (7.7)
qu

 

169

 

\/\.

 

 

 

 

 

 

 

 

 

Figure 7.2: Top panel is the ground state structure of system 1 (without corrugation
potential) and the bottom panel is domain and domain wall structure which is the
ground state structure of system III (with corrugation potential).

170

1.25

1.00

0.75

0.50

0.25

 

0.00 "
0.0 0.2 0.4 0.6 0.8 1.0

Figure 7.3: Dispersion curve for the system I (dotted lines) and the system 111 (solid
lines). Notice the difference not only in the shape but also in the acoustic branches.
The direction of q is along the (10) reciprocal lattice vector of the unit cell.

 

171

where Vj’s are the 2N normal mode frequencies. To get the PDOS, we have carried
out the summation over a regular grid of q values by taking a sufficiently large area
in the q space. If enough points are selected (about 100,000), 9}](V) converges to
the correct result [17]. In Figure 7.4 the result for the system I (without corrugation
potential) is given and provides a check of our calculation procedure because the
ground state of this system is known to be a triangular lattice whose PDOS is exactly
known. The PDOS for system III shows a gap in the low frequency region where
vibrational modes are forbidden due to the commensurate nature of the ground state.
The gap depends on the strength of the corrugation potential. Spiky structure in
gg(V) for system III indicates that its normal mode dispersions are much smaller

than that for the corrugation-free system.

7 .3 Phonon Density of States through Velocity
Auto-correlation Function

For system 11 whose ground state is an incommensurate structure, it is difficult to
calculate the PDOS associated with harmonic vibrations gg(V) by going through
the normal mode frequency analysis. Although the systems II and III differ only
slightly in density, this small difference leads to completely different ground state
structures and it is very difficult to handle the dynamics of the system 11 analytically
because we cannot choose a simple unit cell with a small number of basis atoms for
this system. Strictly speaking, there is no unit cell for this system. The way one can
handle this problem is through MD simulation. It is well known that Z (V) (Fourier
transform of VAF) is directly related to the PDOS in a harmonic system [18]. In
Figure 7.5, we compare Z (V) for system 11 and system III and can see the effects of

corrugation potential and atomic defects on the density of vibrational states. The

172

 

 

 

 

 

 

.. I I fl l I I I I I I I I I l I I Ifi .1
I- d
8 — —-I
u If
- J
- '1
6 b —
A C I
§

V .- .4
.5: ~ -
4 f- '-
.. d

.. I
. i. I 1

I
2 r I“ l —
. d
l- l \ ' d
I \ I
r- , K.~’ d
4'

t’l 1 1 n 14 1 I 1 l 1 1 LL4++u

 

0
0.0 0.5 1.0 1.5 2.0
V

Figure 7.4: Phonon density of states for the same density of system III but without
corrugation potential (dotted line) and system III (solid line) through lattice dy-
namics in harmonic approximation with next nearest neighbor interaction. Notice
the gap for the system III due to corrugation potential.

 

173

 

 

-.I I I I I I I I 1 IT I I I I Iii 1"de
3— 2 ex = 0.9 —
r- :i t d
I- 2 cf
2 z" —

Z(V)

 

 

 

 

lllllJl41lllLlllllllll

O 1 2 3 4
V (THz)

 

Figure 7.5: Fourier transform of the velocity auto-correlation functions near 150K.
Notice the small difference in Z (1!) between system II (solid line) and III (dotted
line) in low frequency region. In system II which is an incommensurate solid, there
are long wavelength phonons which appears to be a linear density of states. Low
frequency modes are remarkably suppressed due to the corrugation potential.

 

174

PDOS for the system II shows the existence of a few allowed modes in the low
frequency region. Physically, the ground state of II can be thought of as that of III
(periodic domain-wall structure) with a few atomic defects caused by the additional
intercalants in the former. We attribute the small difference in the low frequency
PDOS between systems II and III to the existence of long wave length phonons in
the incommensurate solid. The reason we do not attribute this difference to thermal
effects is because the temperature is nearly the same in both the systems. Also, the
number of thermally excited topological defects such as disclination pairs are nearly
the same in both the systems (see chapter 5). The very small tail in Z (V) seen in
the system III in the low frequency region (see Figure 7.5) at this temperature is
due to anharmonicity and/or thermally excited disclination pairs, and major source
of the difference between the systems 11 and III is primarily due to the defects in

the former.

Parameters for the system II were chosen as appropriate for the RbCu. There-
fore, we can compare the results for this system with experimental data. In F ig-
ure 7.5, we see that Z(V) for the system 11 shows a peak (actually a doublet struc-
ture separated by about 0.3THz) near 1THz and its overall width is about 0.5THz.
Furthermore, Z (1!) does not show any discernible peak near 2THz contrary to the
experimental g(V) (see Figure 7.1) given in Ref. [7] obtained from the incoherent
dynamic structure factor. To unravel this puzzle we must look at the procedure KZ

used to obtain g(V) from their inelastic neutron scattering experiment.

 

175

7 .4 Phonon Density of States through Incoher-
ent Dynamic Structure Factor

To understand why KZ saw two pealrs one near 1THz and the other near 2THz,
we briefly review the method K2 [7] used to obtain g(u) from their neutron scat-
tering data. Although Rb is a coherent scatterer, if one takes a sufficiently large q
and considers all possible crystalline orientations in the intercalant plane (perpen-
dicular to the graphite c-axis) then one can use the incoherent approximation i.e.
Swh(q,u) 2 S.(q,u). Basically one is summing over all phonon states in the first
Brillouin zone. The incoherent dynamic structure factor, S.(q, u), can be expanded

as follows [19]:

5.01, u) = 25501.1!)

 

= Z e”W(9) ( In,2 )1" fi U” dugg—gl-1§:—:)-n(ug)) 6(u + 29),), (7.8)

P pl 47"" i=1 ‘°° 6:!
where p is the p-phonon contribution and W(q) and n(u.-) are Debye-Waller factor
and Bose-Einstein factor, respectively. The one- and two-phonon contributions in

Eq. (7.8) are respectively

 

 

_ h 2 n 11 +1
Sim) = e ”(”83", ( ) gum. (7.9)
2 _ -2W 1 I"!2 2
5.01, V) — e ”’5'! (W
°° 9H(V1)9H(-V — V1)
[.00 dlll V1(-V _ V1) n(u1)n(—V — V1). (7.10)

We define a new function, A(q, :1), given by

A(q, u) = «was. ’51—)? (7.11)

176

where c(q) is a normalization constant. This A(q, u) is what KZ referred to as g(V).
If the one phonon term 5'} (q, u) in Eq. (7.8) is dominant, then one can get the PDOS
gH(u) from S.(q, V) i.e.,

n(u) + 1

 

M = A(q. u) z c(q')s:(q. u) = gm). (7.12)

Therefore, A(q, u) or g(u) is the PDOS, 93(0), associated with one phonon excita-
tions when (l) the incoherent approximation is valid and (2) multi-phonon contri-
butions to S,(q, u) are negligible. Thus g(u) is equal to gg(u) only when the above
two approximations are valid. The fact that experimental g(u) and our calculation

930/) = 2Z (V) do not agree suggests that either (1) or (2) (or both) break down.

Now let us study the g(u) for the stage-2 R5024 via MD simulation. We can
calculate directly .S'.(q, V) related to the incoherent scattering cross section; unlike
the experimental measurement we are not restricted to only large values of q to pick
up the incoherent part of the dynamic structure factor. But to compare with the
experiment, we choose q = 4121-1 which is in the range of the q values KZ used.
Figure 7.6 gives A(q,u) (see Eq. (7.11)) for q = 421-1 for may. in both the solid
and the liquid obtained following the procedure of KZ but by using S.(q, u) from our
MD results. We see two peaks with the high frequency peak near 2THz. Our result
for A(q, u) agrees with what KZ plotted as g(u) (circles in Figure 7.6), particularly
very well for the liquid state. (The large intensity seen in the experiment near
3THz is most likely due to the intercalant dynamics perpendicular to the graphite
plane. Also the experimental resolution near 3THz is about 1THz, which implies
large error bars.) In the solid, the doublet structure we obtain near 1THz is washed
out by experimental resolution of ~0.3THz near 1THz‘[7]. The agreement with
the experiment is quite remarkable and tells us two things. One, the incoherent

approximation made by KZ [7] is reasonable for q ~ 4.21"1 and two, our potential

177

 

4IIIIIIIIIIIII‘IITIIIIII.
‘
d

- o L1qu1d :
O Q o
'.... o

TIII

.
O

O. o om

o Cb“. O 00

r, O

 

N
(r'UIIII'fi1"II
llllllllllll

A(q.V)

OHNUIFOIO

0 Solid

II D.(. ”F, &
O o " o

llllllllfillLlIl

0 1 2 3
V (THz)

 

O
l-Jllllllllllllllllllllllll

al.

Figure 7.6: A(q, V) at q = 4A”. The symbols are experimental data [7] for Rng4.
The top and the bottom panels correspond to the liquid and the solid states of
the intercalants respectively. MD simulation was done at T=l44K for the solid
and T=250K for the liquid; the melting temperature in MD simulation is ~ 220K.
Experimental data is for T=100K=0.6Tm for the solid and T_=210K> Tm for the
liquid where experimentally obtained Tm is 162K.

178

is adequate to represent the stage-2 Rb GIC. However, we still have a puzzle, the
PDOS obtained through Z (I!) i.e. 930/) does not show a peak near 2THz in contrast
to the one obtained via S.(q, u) for q = 4A“, the latter of course agreeing with the
results of KZ [7]. To solve this puzzle let us look at the multiphonon contribution

to S.(q, 11).

Since p—phonon contribution to S,(q, u) is proportional to the q“? [see Eq. (7.8)],
one phonon contribution is dominant at small q [20]. Figure 7.7 gives A(q, u) for
q = 0.41214; A(q, V) does not show any high frequency peaks and is a replica of
the Z (11) shown in Figure 7.5. This answers the puzzle i.e. to get the one-phonon
density of states 930/) via the incoherent dynamic structure factor one must confine
A(q, V) to small values of q for which Eq. (7.12) is valid. However it is impossible
experimentally for RbCu since the Rb ions are strong coherent scatterers and for

small q coherent scattering dominates over the incoherent scattering.

We further analyze the PDOS for RbC'u obtained from A(q, u) with q = 0.421“.
To see the effect of the corrugation potential, we calculate gH(V) for a repulsive
screened Coulomb system with the same density as RbCu but in the absence of
graphite corrugation potential. In Figure 7.7, we give the PDOS for both these
systems, the corrugation free system is given in the inset of the figure. We see
clearly that the low frequency modes are strongly suppressed in RbCu due to the
corrugation potential. One important feature of Z (V) or A(q = 0.4, u) is their narrow
width (~ 0.5THz). This narrow width not only gives a well-defined low frequency
peak seen in experimental g(u) but also leads to a pronounced two-phonon peak
when q is large. As we increase the temperature beyond Tm (~ 220K in our MD),
we find that the strength of the low frequency peak starts to decrease and that of

the high frequency peak increases (see Figure 7.6), as expected for a multi-phonon

 

179

 

 

 

 

 

 

     

 

 

 

i i i
. 3 _— I
4 '— I- :
I A Z I
- 3 2 r :1
I- g P '\
3 — Q0 ,. :t II
A : :J\\ l:
h. I- 1 _I \\l |
3 2 Z. :’
4 7 0 Jinnllrnajljnnnlnnnn
I 0.0 0.5 1.0 1.5 20
I-
1 _ V (THz)
0
LLJ 1L1 l L l l l l l l l l l 1 L1 L4
0 1 2 3 4
v (THz)

Figure 7.7: A(q, u) at q = 0.413"1 for the system II at 144K. This agrees very well
with 2Z(V) in Figure 7.5. Inset is the PDOS, gH(I/), for the system I (without the
corrugation potential) obtained from lattice dynamics calculation.

5

180

response. The overall peak near 1THz is quite well defined even above Tm although
the doublet structure is weakened because of the softening of the TA modes. The
reminiscence of the TA mode in the peak is interpreted as a characteristic feature of
lattice-fluid as KZ pointed out. Even if they did not assign two peaks correctly, the
interpretation of the persistence of the TA mode is reasonable. When we increase
the temperature to about twice the melting temperature, the intensities of the one-
and two-phonon peaks become comparable. This is clearly seen in experiments for
both K 024 and Wu [7] thus confirming our multiphonon description of the high

frequency peak.

To see how the relative spectral weights of multi-phonon responses change with
q, we calculate the A(q, V) for a harmonic oscillator with a fundamental frequency
1THz and at temperature 150K. Fig. 7.8 shows that the multi-phonon contributions
become important with increasing q as expected (Eqs. (7.9) and (7.10)). The domi-
nant peak at 1THz is due to the one-phonon response. We see the evolution of two-
and three-phonon contributions at 2THz and 3T Hz, respectively as we increase q.
This clearly tells us that to get the PDOS we must confine A(q, V) to small values

of q.

Finally, let us calculate multiphonon response of the system I which is a corruga-
tion free system and compare it with that of the system II. Figure 7.9 shows A(q, V)
at q = 413'1 for the systems I and II both in the solid and in the liquid state.
For the system I, either in the solid or in the liquid state, we cannot distinguish
one-phonon from multiphonon contributions in A(q, V). This tells us that a clear
distinction between one and multiphonon process in A(q, V) is indeed a remarkable

characteristic of strong corrugated system (with a narrow PDOS).

 

181

 

A(q.V)

 

 

 

Kg

 

 

 

Figure 7.8: Frequency dependence of A(q, V)’s for different values of q from 0.4 to
4 in steps of 0.4 (in units of A") for a simple harmonic oscillator of fundamen-
tal frequency 1THz. This clearly shows the evolution of multi-phonon peaks as q
increases.

182

 

 

 

 

 

F..." 1 I I I I I I I r T I I 1 I I
”5 E' .. ...,. Liquid
C 3" 's.
0.04 _— :'. .‘t,’
E x
0.03 .— , ‘1.
0.02 E- ’
I; 0.01 :- .
¢ : .. j l l l l I l l I l l l
I' 0.00 p I Ifi I : I I I I f I I I I E I I
if 0.10 _—
: Sofid
0.08 :-
E i
0.06 L- 'y -_
: f ‘
0.04 :- g é
- f '-=.
0.02 E— ‘: “tax
0.00 ". 11 Lik1 1 1L11T‘1“ “1

O

1 2 3
V (THz)

Figure 7.9: A(q,V) at q = 4A“ (multiphonon response) of the systems I (dotted
line) and II (solid line). The system I (system II) corresponds to at T=120K (144K)
for the solid and T=233K (250K) for the liquid; The melting temperature is about
200K and 220K for systems I and II, respectively.

183

7 .5 Summary and Conclusion

In conclusion, we have studied the phonon density of states for the stage-2 Rb GIC
by using MD simulations. The advantage of the MD simulations is that we can
obtain both the incoherent dynamic structure factor and the Fourier transform of
the velocity-auto correlation function independently and by comparing these two, we
can give a clear interpretation of the origin of the two peaks in the density of states
obtained from the incoherent dynamic structure factor at large values of q. The
results of the inelastic neutron scattering experiments of Kamitakahara and Zabel
[7] for this system have been reproduced by MD simulations and we have clearly
demonstrated that the multi-phonon response is extremely important at large q
particularly for systems with strong substrate corrugation. We also predict that
the relative strength of the high frequency peak compared to the low frequency one

should increase with q which can be easily checked experimentally.

Bibliography

[1] R. J. Birgenau and P. M. Horn, Science 232, 329 (1986); F. F. Abraham, Adv.
Phys. 35, 1 (1986).

[2] Graphite Intercalation Compounds, vol. 1, Structure and Dynamics,
ed. by H. Zabel and S.A. Solin, Springer Series on Topics in Current Physics,
Springer (1990).

[3] S. N. Coppersmith et al., Phys. Rev. B 25, 349 (1982).
[4] Edward Vives and Per-Anker Lindgard (preprint).
[5] W. A. Steele (preprint).

[6] M. Nielson, J. P. McTague and L. Passell, Phase transitions in Surface
Films, edited by J. G. Dash and J. Ruvalds, p127 (1979).

[7] W.A. Kamitakahara and H. Zabel, Phys. Rev. B 32, 7817 (1985).
[8] H. Zabel in reference [2], p101.

[9] M.M. Bredov, B.A. Kotov, N.M. Okuneva, V.S. Oskotskii, and A.L.Shakh-
Bugadov, Fiz. Tverd (Lennigrad) 9, 287 (1967)[Sov. Phys. Solid State 9, 214
1967”.

[10] Hyangsuk Seong, Surajit Sen, Tahir Cagin, and S. D. Mahanti, Phys. Rev. B45,

8841 (1992).
184

185

[11] R. Clarke, J. N. Gray, H. Homma, and M. J. Winokur, Phys. Rev. Lett. 47,
1407 (1981); H. Zabel, et al., Phys. Rev. Lett. 57, 2041 (1986).

[12] W.A. Kamitakahara, N. Wada, and S.A. Solin, Solid State Commun. 44, 297
(1982).

[13] J. D. Fan, George Reiter, and S. C. Moss, Phys. Rev. Lett. 64, 188 (1990).

[14] Hyangsuk Seong, S. D. Mahanti, Surajit Sen, and Tahir Cagin, Phys. Rev. B46,
8748 (1992).

[15] Hyangsuk Seong and S. D. Mahanti (unpublished).

[16] Quantum Theory of the Solid State; Part A, Joseph Callaway, Academic Press
(1974).

[17] Garboczi, Ph.D dissertation, Michigan State University (1984).
[18] A. Rahman, M.J. Mandel and J .P.McTague, Jour. Chem. Phy. 64, 1564 (1976).

[19] W. Marshall and S.W. Lovesey, Theory of Thermal Neutron Scattering,
(Oxford University Press, Oxford, 1971).

[20] RA. Egelstaff, Dynamics of Disordered Materials, edited by D. Richter,
A.J. Dianoux, W. Petry, and J. Teixeira, p2 (1988).

Appendix A

Average length for A 32% 0

A.1 6 - bonds (same layer)

(LIMA 34A = (L); +a;83:A
+ 4(535 + 253..)$u(L94.4 - L043)

"' 23ul2xub36 + (xii-1 + “pi-”(bk + bZnHA

(Lum EB = (L); + 025(1 _ 33»le
+ 4(935 + 253..)(1 " 9"»)(11944 - L33)

_ 2(1 - zu)[2$ubis + (“in-1 'l' 3u+1)(b;u + bEnHA

me)?” = (L); + 035241 — WA
+ 4(bz. + 2b;.)(1 - 22..)(L‘3... - Lie)
- 2(l - MHZ-tuba + (me-1 + xu+1)(b;V + b3.)]A
l+m

(Lls = T-l-I—m leul’gA + (1 - xungB — $u(1 - xulAl .

u=1

Q-5
2)

 

a3, = 4%[a B (or1 .3] sin2(

186

(A.1)

(A.2)

(A.3)

(AA)

(A.5)

187

5;, = §[3-B(er~V]sin(Q-3)sin(Q-I‘x) (A.6)
.;,, = glv-er‘-i1(1—cos(Q-a)(1—oos(Q-i) (Arr)

A.2 11(1)) - bonds (interlayer)

(Lam-lltA = (Lam-1).. + 0"(1’u-1 + quLSA "' L33)
+A [ a;u($u-23u-I - III-1% + tutu“)
- 25;»(31-1 + 3,24)
-‘ 25:u(-'€u—2$u-1 + zflzfl‘l'l)
- 53.,(1'u—23u-1 + 2344-le + $u3u+ll

+ 8:.(zu—2ze—1 - helm. + mixed] (Ao8)
(Lam-llfa = (Lum-lly + a"(1 — “in-1 + 1 - 2..)(13‘1. - Lie)
+Ala:.((1- xe-2)(1- zu-n) - (1 - zu-1)(1- xe)+(1- xe)(1 — n+0)
+2b3..((1 - xu-llzu-l + (1 - mare) + 2b:..(a=..-2(1 - ace-x) + (1 - xulxu-H)
+b:.((x..-2 + xu)(1 - xe—1)+(1 - axe.-. + $u+1))
-s:.((1- are-2X1 - xu-I) - 2(1— rep-00 — 3V)+(1— we)(1— x..+1))]
(A.9)

(Lum-llfiB = (Lum-lly + an“ '— xu-l — $u)(L91A - L948)

+ A S (A.10)

(1 "' zulxu-l “I" $u(1 — zu-I)

 

where

188

5 = a:u($u-2$u-I — “Bu—13a + III-tun)

+2534“ " 23a)“ — zu-llzz-l + (1 — 23nd)“ - tulxi.)

+2b:y($u-2$u-I(1 - ire—1X1 " 23») + (1 - Zia-1X1 - xulxuzuH)

+b:n((3u+l + “Bu—1X1 - a=u)=¢=u(1 - 233nd) + (tn-2 + In)“ " 3»-I)$u-I(1 - 2%))
--9:.,((3u+1 + “Bu—1X1 " tal-WU " 22,._1) "' (tn-2 + $00 " zu-llxu-lu - 2%))

(A.11)

Appendix B

EMT in Harmonic Limit

In the harmonic limit where 7.4 = 13 = 0, Eq. (3.32) reduces to a simple form,

F th. + Kali;

 

 

= —— B.l
U“) K; + K... K; + Kg. ( )
Eq. (3.28) then becomes
F
= —. B.2
0») h. + K. ( )
Eq. (3.33) gives
F _ F Kéh, + thg
h‘+ K. " (1—$)K£+KA H1”) K;+KA
F K’h, K ho
+ ° + B 1’1. (3.3)

 

mm + a: K; + KB
This equation is satisfied under an arbitrary force F. If we rearrange Eq.(B.3), we
have an equation such as aF + ,6 = 0. Solutions of this equation for the arbitrary
force F are a = 0 and fl = 0. Therefore, we obtain two equations which are exactly

the same as the results of Thorpe et al. [14].

189

Appendix C

Mean Square Displacement in
d-dimension

Let me show here how Eq. (6.6) is obtained. All what we need to know is a volume
element dr in d—dimensional spherical coordinatel. To get the dr, we use a following

identity.

00 dam"1B2 0° dye"? . - - 0° dze"a = [00° drSdrd'le”: (C.1)

where dr 5 Syd-1dr and r2 in d—dimension is a:2 + y2 + . ~ - 22.
L.H.S. of Eq. (C.1) = 1rd”.
R.H.S. of Eq. (C.1) =SdI‘(d/2)/2

where F is the gamma function. Therefore, 54 = 2fid/I‘(d/2).

Now, we are ready to calculate the mean square distance in d—dimension.

 

 

<(r(t>-r(o»’> = fdrr’a.(r,t) (0.2)
4/2
d/2 '
= (aim) S‘Idrrwe—rmam (0'4)

 

1Shang—Keng Ma, Statistical Mechanics, World Scientific, p77 (1985) ; M.F. Thorpe, Lecture
note

190

191

 

= (217275)”: 54 / fig—his 2a(t) “ed/2e” (C.5)
= (wall/25 «Tux/20m“ I‘M/2+1) (0.6)

a(t)d (C.7)