I}!
’£:I!.-a..l.2i.r
Siélufl
p 2. $

.14..
I: . . .

(x

.f it: v’.’;.ll§‘ru.hlua.iiflvart..t
I? if! 1....
42?. 3.112).}!

(0...?

..¢.I‘u§v‘...lutir.v
(agar-livid]: {glitz
11$»..lfan1lxll4lly‘

£373:

5.2.33.1..11
‘ 1.3.1. (2.2L! .3. .a
.1

1552...! 1. v . .‘
a.:l..!..l1ti€lll!!
~zIAI.J.l...l:...y4/I: .31:
((113 :3 la 1

 

THESIS

TA E UNI

TI lung

93()

lllll'llllllllllllll

904 5877

 

 

 

will

This is to certify that the

dissertation entitled
CURVATURE ENERGY OF NUCLEI,
ELECTRONIC EXCITATIONS OF CARBON CLUSTERS,

AND THERMAL PROPERTIES OF SODIUM CLUSTERS
presented by

Nengjiu Ju

has been accepted towards fulfillment
of the requirements for

Ph.D. Physics

degree in

 

 

8850

Date July 20, 1993

 

MSU is an Affirmative Action/Equal Opportunity Institution 0-12771

 

LIBRARY
MIchIgan State
University

 

 

 

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

DATE DUE DATE DUE DATE DUE

 

 

 

3
H1213

'E

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU Is An Affirmative Action/Equal Opportunity Institution
rc\daledue pmlp. I

CURVATURE ENERGY OF N UCLEI,
ELECTRONIC EXCITATIONS OF CARBON
CLUSTERS, AND THERMAL PROPERTIES

OF SODIUM CLUSTERS

By

Nengjiu J u
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

CURVATURE ENERGY OF NUCLEI, ELECTRONIC
EXCITATIONS OF CARBON CLUSTERS, AND THERMAL
PROPERTIES OF SODIUM CLUSTERS

BY

Nengjiu Ju

This thesis consists of three separate parts. In the first part we calculated the
nuclear correlation energy EMT in the RPA approximation with particular attention
paid to the A-dependence. We found that Eco" has a very different A-dependent
behavior from the mean-field total ground state energy Ema”. In this way, we want
to shed some new light on the longstanding contradiction between the curvature
coefficient value ac resulting from mean-field calculation which is always close to 10
MeV and that from the direct adjustment of mass-formula to nuclei ground state

masses which is zero.

Part two is dedicated to the study of the multipole electronic collective excitations
in carbon clusters using linear response theory. The single—particle wave and energy
levels are determined by a tight—binding Hamiltonian, whose parameters have been
determined by a fit to the LDA results. Because 060, the so—called buckyball, has a
spherical symmetry, its treatment is much simpler than that of an arbitrary cluster,
it is naturally our first choice to apply our model. With some work, we were also
able to completely solve the RPA problem within the tight—binding framework for an

arbitrary geometry of the cluster. We applied this new formalism to other carbon

clusters. The equivalents of both 71' and a plasmons in graphite are predicted for
these clusters. Since only dipole excitation is optically allowed, the other multipole
excitations are best studied by electron energy—loss spectroscopy (EELS). To this
end, we calculated the differential cross sections of inelastic electron scattering off
these carbon clusters. We compared our dipole electronic excitations with available
photoionization experimental data and our differential cross sections with available
EELS data. Good agreement is obtained between our calculations and available

experimental results where they are available.

The third part addresses the thermal properties of sodium clusters using the
isothermal molecular dynamics developed recently by Aurel Bulgac and Dimitri Kus-
nezov. A wide range of properties of small sodium clusters are investigated and pre-
sented. We found that sodium clusters undergo two “phase transitions”, one around
200 K from a crystal to a glassy or molten state, and a second one around 800 K,
from a molten to a fluid state. At low temperatures these clusters are essentially
incompressible, but relatively easy to deform. At high temperatures they become

extremely soft and evaporation of atoms sets in.

To my wife, Wenge Zhang, for her love, understanding, tolerance and to our

daughter, Julie Ann Ju, for her laughter.

iii

ACKNOWLEDGMENTS

My deepest thanks go to my adviser, mentor and friend, Professor Aurel Bulgac,
for his patience and encouragement during the past five years, for his kindness and
humanity, for his understanding and supportive altitude. He has taught me not only
physics, but the values of hard work and devotion. In a word, I am very grateful to

him.

Thanks also go to my guidance committee members, Professor W. Benenson,

Professor G. F. Bertsch, Professor S. D. Mahanti and Professor W. W. Repko.

Thanks are also due to members of our weekly cluster meeting: Professor D.

Tomanek, Dr. Yang Wang, Dr. C. Lewenkopf and Mr. V. Mickrjukov.
Thanks are also due to all my friends here.

I wish to express my special thanks to members of our host family: Margaret
Houseman and David Houseman, and their two lovely children Loran and Sara Maria,

for their providing us numerous happy and memorable gatherings.

iv

Contents

LIST OF TABLES
LIST OF FIGURES

1 Introduction
1.1 Nuclear Correlation Energy .......................
1.2 Atomic Clusters ..............................
1.2.1 Multipolar Responses of Fullerenes ...............
1.2.2 Thermal Properties of Sodium Clusters .............
1.3 Thesis Organization ............................

2 Correlation Energy and Curvature Energy of Nuclei
2.1 Introduction ................................
2.2 RPA Formalism of Correlation Energy .................
2.3 Numerical Details and Results ......................

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

3 Collective Electronic Excitations in Fullerenes
3.1 Introduction ................................

3.2 Theoretical Description ..........................

3.2.1 Description of the Tight—binding Hamiltonian .........
3.2.2 Linear Response Theory of the Tight—binding Hamiltonian . .
3.2.3 Outlines of the Calculations of C60 ...............

3.2.4 Outlines of the Calculations of C20,70,100 ............
3.3 Theoretical Results and Comparison with Experimental Data .....

3.4 Discussions ................................

vii

viii

3O
30
31
32
34

37
39
54

— I! ”an ~ inn», '-— - .
—-A-.._‘..«-.-.. ._ _‘_. - my... . .- ~.

 

4 Finite Temperature Properties of Small Sodium Clusters
4.1 Introductory Remarks ..........................
4.2 Isothermal Dynamics ...........................
4.3 Thermal Properties ............................
4.3.1 Thermodynamic Properties ...................
4.3.2 Geometric Properties .......................

4.4 Closing Remarks .............................

5 Summary and Perspectives
5.1 Summary .................................

5.2 Perspectives ................................
A Multipolar Responses for C60

B Multipolar Responses for C20,70,1oo

vi

102

108
108
109

111

118

List of Tables

2.1 Correlation energies for magic number nuclei .............. 27
2.2 Mass—formula fit of the magic nuclei correlation energies ....... 27
2.3 Correlation energies for open—shell nuclei ................ 28
2.4 Mass—formula fit of the magic nuclei correlation energies ....... 28
3.1 Dipole oscillator strength of C60 ..................... 49

 

List of Figures

1.1
1.2
1.3
1.4

3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
3.10
3.11
3.12
3.13
3.14
3.15

4.1
4.2
4.3
4.4
4.5
4.6

Mass abundance spectrum of NaN clusters ...............
Structure of the C60 “buckyball” cluster ................

Giant dipole resonance of C60 ......................

photoionization yield of Cgo ....................... 10
Tight—binding spectrum of the C60 ................... 33
Multipolar responses of C20 ....................... 40
Multipolar responses of C60 ....................... 41
Multipolar responses of C70 ....................... 42
Multipolar responses of C100 ....................... 43
Ionization spectrum of C60 ........................ 44
Ionization spectrum of C70 ........................ 45
Dipole response of C100 along the symmetry axis ............ 47
Dipole response of C100 perpendicular to the symmetry axis ..... 48
The dynamical structure factor of C60 .................. 52
The dynamical structure factor of C70 .................. 53
Electron energy loss spectroscopy of C60 ................ 55
Electron energy loss spectroscopy of C70 ................ 56
Feynman diagrams of the a —-> 2 7r and a -—> 3 7r processes ....... 59
Feynman diagrams of the decay of a plasmon through the emission of

2 or 3 electrons .............................. 60
Density of states of N as, 14, 20, 30, 40 .................... 80
Potential energy as a function of temperature ............. 84
Specfic heat as a function of temperature ................ 85
Relative bond length 6 as a function of temperature .......... 87
Pair correlation function for Nago .................... 89
Pair correlation function for N am .................... 9O

viii

4.7 Allowed region of fl and ”y ........................ 92
4.8 Principal momenta of inertia ....................... 93
4.9 Relative covariances of principal momenta of inertia .......... .94
4.10 RMS radius and its covariance ...................... 96
4.11 Shape parameter [3 and its covariance .................. 97
4.12 Shape parameter 7 and its covariance .................. 98
4.13 Ion density distribution for Nazo ..................... 100
4.14 Ion density distribution for Na“) ..................... 101

 

Chapter 1

Introduction

1.1 Nuclear Correlation Energy

One of the well—established approaches to obtain nuclear binding energies is based
on the analogy of a nucleus to a drop of incompressible fluid. The first evidence in
support of such a simple “liquid drop” model for the nucleus comes from the fact
that per nucleon binding energy is about the same for all nuclei. If we are willing
to ignore the small local departures, it is possible to develop a simple formula that
expresses the binding energy of a nucleus in terms of its number of nucleons (in other
words, its volume since the volume is proportional to its number of nucleons) along

with surface correction, curvature correction and higher—order corrections.

The liquid drop model (LDM) has been very successful in determining many nu—
clear properties. An expansion of the nuclear energy of a finite nucleus into volume,
surface, curvature and higher—order contributions has proved enormously useful in
calculating fission barriers, ground—state masses and deformation, and other nuclear
properties [1, 2, 3]. However, a serious anomaly exists between the value of the
curvature energy coefficient aC obtained theoretically and that from adjustments to

experimental data [4]. Direct adjustment of mass—formula to the ground state masses

of nuclei throughout the periodic table yields values of the curvature energy coeffi—
cients aC (aC in the term aCA1/3 in the mass-formula) close to zero. On the other
hand, many previous theoretical calculations have given a typical value 10 M eV for

do.

By summing up all the ring diagrams, we calculated the nuclear correlation energy
ECO“. in the RPA approximation with particular attention paid to the A—dependence.
We found that Eco" has a very different A—dependent behavior from that of the
mean—field total ground state energy Emean. In this way, we want to shed some new

light on the anomaly aforementioned.

1 .2 Atomic Clusters

Clusters are aggregates composed of a countable number of atoms, which can be
as small as two or as large as hundreds of thousand. The real breakthrough comes
through the discovery of quantal shell structure in small droplets of sodium metal,
with characteristic “magic numbers”, by Knight and coworkers [5] , see fig. 1.1.
Since then, cluster physics has been an active field of theoretical and experimental
investigation. Many properties of a large variety of clusters have been examined
[6] — mass abundance spectra, fragmentation spectra and binding energy, supershell
structure [7], ionization potential [8, 9], photoelectron spectra and electron affinity,
static [11, 12, 13] and dynamical [14, 15, 16] polarizability, plasmon resonance spec-
tra [17] and thermal properties [18, 19], to name just a few. With the advent of
new technology, it has become possible to have well—controlled cluster sources, thus
to make clusters of well-defined sizes and to measure their properties precisely. Be-
cause clusters lie somewhere between a solid and an atom, they provide for a test

of our understanding of bonding and of the transition from atomic to molecular and

 

 

 

 

 

 

 

 

 

 

 

 

:3
(a)
'.‘.' '.‘.'.'.'}.0.‘.‘.‘.-'.l.'.l.’.'.'m-'.‘.‘
2 L
3 92
g
I] [l
L L l 1 r '
1.6 m
.19 (b) O 4 1.
1.2 - '
3 .
E. 0.3 _
4
A D
f: 0.4 -
a I
4 0 _
-0.4 WM
3 20 34 40 5' m
m ct scan m not". N

Figure 1.1: Mass abundance spectrum of NaN clusters. (a)Mass abundance spectrum
of NaN clusters, N:4-75. (b) The electronic energy difference. A(N + l) — A(N) vs
N. The labels of the peaks correspond to the closed—shell orbitals. [From Knight et
al, Phys. Rev. Lett. 52, 2141 (1984)].

 

 

 

 

ultimately bulk properties. A compelling stimulus has been the desire to understand
how an extended crystalline solid develops from growing cluster aggregates. Clusters
possess many unique properties and they can be used to make new materials. Stimu—
lated perhaps by its potential use, from its very beginning, cluster physics has been a
very active field, bringing methods and concepts from atomic and condensed matter
physics, quantum chemistry, thermodynamics, surface science, optics, nuclear physics

and others into unprecedented combinations.

The most thoroughly studied clusters are the sodium clusters. Since sodium clus-
ters can be generated in conventional sources under supersonic adiabatic expansion
conditions, they can be studied for a longer time than it is possible for most other
clusters that require the use of the Smalley laser vaporization source [43]. The first
evidence of shell structures in clusters comes from the study of sodium clusters [5].
Recently, abundance spectrum has been measured up to 22,000 atoms and has re—
vealed interesting new shell structures [45]. The persistence of the strong shell effects
up to very large numbers of atoms in a cluster is a quite remarkable phenomenon

[46, 45].

The most interesting clusters studied so far are certainly the carbon clusters,

especially the so—called fullerenes including the celebrated buckyball C60.

To explain the large abundance of clusters with 60 carbon atoms, the existence of
a stable, spherical C60 molecule, the so—called buckyball, has been proposed [44]. In
the suggested structure, twelve regular pentagons and twenty hexagons are connected

to form an icosahedral structure, see fig. 1.2.

The most exciting development of the study of C60 has been the discovery of the
fourth form of pure carbon solid, the so—called fullerite, which was discovered about

only three years ago if we consider the amorphous carbon also as a pure solid form

 

3"“. ’:; :5‘

 

 

Figure 1.2: Structure of the C60 “buckyball” cluster.

(the other two forms are sp3—bonded diamond and spz-bonded graphite). In the fall
of 1990, W. Kratschmer et al [47] made a major breakthrough in the synthesis and
separation of C60. By chemical extraction methods, they obtained large quantities of
C60. This breakthrough has led to the intense experimental and theoretical studies
of the structural [48, 49, 50, 51] and electronic properties [52] of this novel solid form

of carbon.

The most striking property of fullerite is perhaps its superconductivity when it is
doped. Although the highest critical temperature of doped fullerite is still well below
that of the the high—TC superconductors, these intercalated fullerites have polarized
the interest of many investigators because they have superior material properties and

hence bear a higher potential for applications.

Many theories have been developed to explain a large variety properties of different
clusters, for example, the magic numbers of alkali metals. The simplest of these is
the spherical jellium model [6] of alkali metal clusters. Because of the delocalized
character of the valence s—electrons of alkali atoms, they are shared by all ions in a
cluster and this feature leads to shell effects similar to those observed in atomic nuclei.
To a very good approximation, the ions can be thought of as simply providing an
almost uniform positive background, and a jellium approximation is quite reasonable
[69]. The valence s—electrons move in a common mean field and the ensuing strong
shell effects account very well for both the enhanced stability of the magic clusters,
shapes corresponding to different numbers of atoms and the character of the optical
response (Mie resonance), which has analogous features to the well known giant dipole
resonance in nuclei. Spheroidal or ellipsoidal extensions of the jellium model are
derived in order to study deviations from sphericity and are used to provide general

shapes of cluster as oblate or prolate ellipsoids [22]. The spheroidal jellium model

also shows that within the jelluim model, special stability at shell closing is associated

with vanishing distortion from sphericity.

Among the more sophisticated methods are the density functional theory (DFT)
based on the local density approximation (LDA), configuration interaction (CI) cal-
culation, and calculation based on the Car—Parrinello method, to just name a few. In
these studies, it has been found that generally each cluster has several isomers at T
= 0 K with a very similar total energy, therefore thermal effects may have important

implications about the properties of the clusters [73].

Although there are many theoretical studies on a variety of different clusters, we
will concentrate on two studies of two classes of clusters in this thesis. The first
one is the study of the electronic excitations in fullerenes, which is presented in
chapter 3. Fullerenes are a group of carbon clusters with hollow cage structures, with
only pentagons and hexagons on their surfaces. The second one is the study of the
thermal properties of small sodium clusters using the newly developed isothermal
molecular dynamics [18, 19, 40, 39], which is the subject of chapter 4. Although
many theoretical studies have been done on clusters, most of them are concerned with
ground state properties, such as ground state geometrical and electronic structures.
Without justification, these calculations have been compared with experimental data.
On the other hand, clusters are produced at high temperatures [41]. Therefore the
study thermal properties of clusters is very important. The choice of sodium clusters

is pragmatic. Sodium clusters have been most thoroughly investigated [5].

1.2.1 Multipolar Responses of Fullerenes

Fullerenes are carbon clusters with hollow cage structure, with only pentagons and

hexagons on their surfaces. The best known member of fullerenes is the aforemen-

tioned buckyball C60. Since fullerenes are essentially hollow graphite shells, they
should support multipolar plasmon oscillations, closely related to the plasmons in
graphite. G. Barton and C. Eberlein [57] estimated the frequencies of these multi-
polar plasmons by adapting the reasonably successful hydrodynamic approach famil-
iar from bulk and surface plasmons [58, 59], with just one parameter calibrated on
graphite. They predicted 7r plasmons in the range between 6 and 8 eV, and a plas-
mons near and above 25 eV. In a more serious and physical treatment of the dipole
plasmon of C60, Bertsch et al [60] used the linear response theory (RPA) to study the
collective plasmon excitations in C60 clusters, they found that the valence electrons
are quite delocalized and show collective excitations. They predicted a giant dipole
resonance at an unusually high energy of 20 eV, see fig. 1.3. Their prediction is some-
what lower than the prediction 25 eV from the hydrodynamic model. The existence
of a giant dipole resonance in C60 was subsequently confirmed by an independent
photoionization experiment on a C60 gas target [61], see fig. 1.4. In computing the
matrix elements of the residual Coulomb interaction arising from the I = 0 transitions
in Ref. [60] for the dipole response, it was assumed that the spatial extent of the s—
and p-states is vanishingly small, this approximation results in a slight overestimate
of the Coulomb matrix elements. Another approximation in ref. [60] was to neglect
the l = 2 component in the particle—hole state [1912 > resulting from a p—hole state to
a p—particle state transition. These approximations make their approach unsuitable
to study other multipolar responses and also the multipolar responses of other non-
spherical fullerenes. It is part of our interest in chapter 3 to extend their calculation
without these approximations and to see to what extent their approximations are

valid. It is also our interest to extend the calculation to other multipoles.

In chapter 3, besides studying the multipolar responses of the C60, we also study

 

BO

 

80 r . . .
70 r "'.........: 7o
60

SO

30'
20
10

8
Infiuans 4019111080 Pa]uJfla]u‘

 

 

 

0 510152025303540
8(0V)

Figure 1.3: Dipole response of C50 clusters to an external electromagnetic field [From
G. F. Bertsch, A. Bulgac, D. Tomanek, and Y.Wang, Phys. Rev. Lett. 67, 2690

(1991)].

 

 

 

 

 

 

   

 

     

m
1.0
0.8 < 40
0.6 '4
04- 35 m.
0.2 . 30 E
0.0 - :1
. . . 25 '
70 7.4 \ 7.8 p ‘5
20« \ 20 ‘°
\ \\
10« 10 2
Q)
on)
5~ s E
0- ' ' Y r ' I r r mm r0
0 4 12 16 20 24 28 32 36

 

photon energy / eV

Figure 1.4: Observed C3}, photoion yield as a function of photon energy. [From
I. V. Hertel et al, Phys. Rev. Lett. 68, 784 (1992)].

11

the multipolar responses of other fullerenes. We have two aims. First of all. we are
interested in seeing how the plasmons vary with the sizes of the clusters. In the
classical picture of Mie, the energy of the plasmon is independent of the size of the
cluster and is totally determined by the free electron charge density of the cluster.
Secondly, we want to know how the plasmons behave with nonsphericity. Since other
fullerenes do not have spherical symmetry as C60 does, each plasmon of C60 should
split in other fullerenes. The extension from the RPA calculations of C60 to those of
other fullerenes is not trivial. Usually the nonsphericity of the system under study
presents a formidable technical problem in RPA, fortunately we are able to solve
the RPA equation for any geometry of the system considered in the tight—binding

framework for the single particle description of the electrons

Since only dipole excitation is optically allowed, the other multipolar excitations
are best studied by electron energy—loss spectroscopy (EELS). For this reason, we cal-
culate the differential cross sections of inelastic electron scattering off the fullerenes.
We compare our dipole electronic excitations with available photoionization experi-
mental data and our differential cross sections with available EELS data. Good agree-
ment are found between our calculations and available experimental results where they

are available.

1.2.2 Thermal Properties of Sodium Clusters

There are several reasons why it is interesting to study finite temperature properties
of atomic clusters, in particular of simple metal clusters. The most obvious one is the
fact that experimental results obtained so far are for relatively hot clusters [6, 74].
Even though the experimentalists do not have at the moment a way of determining

the temperature of the clusters they produce, and only relatively rough estimates are

12

available, it is very likely that these clusters are melted. An estimated temperature
of T m 500—600 K seems to be realistic under experimental conditions [41]. On the
other hand, most of the theoretical studies are performed at zero temperature. With-
out further justification, these results have been compared with experimental data.
Therefore properties and the behaviour of finite systems at nonzero temperatures are

of unquestionable theoretical interest.

We will try to elucidate to some extent the role of a finite temperature in sodium
cluster in chapter 4. Although it is unquestionable that the most ‘fundamental”
way to describe atomic clusters at zero as well as at finite temperatures is ab initio
calculations, the computer time required to extract relevant physical information is
appalling. It is for this reason that we will adopt a less ambitious approach. We
will treat the electrons implicitly by using an effective many—body—alloy (MBA )
Hamiltonian. The parameters of this Hamiltonian were determined from ab initio

calculations with no free parameters [75].

1 .3 Thesis Organization

In chapter 2, we use RPA to calculate the nuclear correlation energy to address the
longstanding contradiction between the curvature coefficient value aC resulting from
mean-field calculation and that from the direct adjustment of mass—formula to nuclei

ground state masses.

We discuss the multipolar responses and differential cross sections of inelastic
electron scattering of C60 and other fullerenes in chapter 3. Comparison is made
between our calculations and available experimental data and good agreement is

obtained.

We elucidate the thermal properties of sodium clusters in charter 4. The electrons

13

are treated implicitly by a MBA Hamiltonian. The possible link between the red shift

of the Mie plasmon in sodium clusters and their thermal expansion is discussed.

A summary of the results obtained and how to extend the works discussed in this

thesis are presented in chapter 5.

Finally the appendices are an integral part of this thesis, they provide the for- _
mulas used. We discuss the formulas of the multipolar responses and EESL of C60
in appendix A and their extension to nonspherical fullerenes in appendix B. This
extension is nontrivial and the formulas are much more complicated than those in

appendix A.

 

Bibliography

[1] W. D. Myers, and W. J. Swiatecki, Ann. of Phys. 55 395 (1969) .
[2] Y. U. Chu, B. K. Jennings, and M. Brack, Phys. Lett. B 68,407 (1977).
[3] M. Brack, C. Guet, and H. B. Hakansson, Phys. Reports 123, 275 (1985).

[4] W. Stocker, J. Bartel, J. R. Nix, and A. J. Sierk, Nucl. Phys. A 489, 252
(1988).

[5] W. D. Knight, K. Clemenger, W. A. de Heer, W. A. Saunders, M. Y. Chou,

and M. L. Cohen, Phys. Rev. Lett. 52. 2141 (1984).

[6] W. de Heer, W. Knight, M. Y. Chou, and M. L. Cohen, in Solid State
Physics, Vol. 40, edited by H. Ehrenreich and D. Turnbull (Academic, New
York, 1987).

[7] H. Nishioka, K. Hansen, and B. Mottelson, Phys. Rev. B 42, 9377 (1990).
[8] D. Snider, and R. Sorbello, Solid State Commun. 47, 845 (1983).
[9] D. Beck, Solid State Commun. 49, 831 (1985).

[10] A. Rubio, L. Balbas, L. Serra, and M. Barranco, Phys. Rev. B 42, 10950
(1990).

[11] D. Snider, and R. Sorbello, Phys. Rev. B 28, 5702 (1983).
14

[121
[13]
[141

[151

1161

v——‘
,_.a
_]

__1

1181

1191
1201
1211
1221
1231
1241
1251

[261

 

15

D. Beck, Phys. Rev. B 30, 6935 (1984).
V. Kresin, Phys. Rev. B 39, 3042 (1988).
W. Ekardt, Phys. Rev. Lett. 52, 1925 (1984); Phys. Rev. B 32. 1961 (1985).

C. Yannouleas, R. Broglia, M. Brack, and P. Bortignon, Phys. Rev. Lett.
63, 255 (1989).

v. Kresin, Phys Rev. B 40. 12507 (1989).

K. Selby, M. Vollmer, J. Masui, V. Kresin, W. A. De Heer, and W. D.

Knight, Phys. Rev. B 40, 5417 (1989).

A. Bulgac and D. Kusnezov, Phys. Rev. Lett. 68, 1335 (1992); Phys. Rev.

B 45, 1988 (1992).

N. Ju and A. Bulgac, Phys. Rev. B July 15, (1993).

O. Sugilo, and H. Kamimura, Phys. Rev. Lett. 65, 2696 (1990).

S. W. Wang, L. M. Falicov, and A. W. Searcy, Surf. Sci. 143, 609 (1984).
K. Clemenger, Phys. Rev. B 32, 1359 (1985).

R. Kamai, and J. H. Weare, Phys. Rev. Lett. 65, 80 (1990)

J. P. Bucher, P. Xia, and L. A. Bloomfield, Phys. Rev. B 42, 10858 (1990).
E. Engel, J. P. Perdew, Phys. Rev. B 43, 1331 (1991).

V. Bonaéié—Koutecky, P. Fantucci, J. Koutecky, Chem. Phys. Lett. 146, 518

(1988); ibid 166, 532 (1990).

 

16

[27] I. Boustani, W. Pewestorf, P. Fantucci, V. Bonaéié-Koutecky, and J.

Koutecky, Phys. Rev. B 35, 9437 (1987).
[28] I. Boustani, J. Koutecky, J. Chem. Phys. 88, 5652 (1988).

[29] V. Bonaéié-Koutecky, I. Boustani, M. Guest, J. Koutecky, J. Chem. Phys.

89, 4861 (1988).

[30] V. Bonaéié-Koutecky, P. Fantucci, J. Koutecky, Phys. Rev. B 37, 4369

(1988).

[31] J. Flad, G. Igel, M. Dolg, H. Stoll, and H. Preuss, Chem. Phys. 75, 331

(1983).
[32] J. L. Martin, R. Car, and J. Buttet, J. Chem. Phys. 78, 5646 (1983).
[33] J. L. Martin, J. Buttet, and R. Car, Phys. Rev. B 31, 1804 (1985).
[.34] J. L. Martin, 2. Phys. D 12, 347 (1989).
[35] R. Car, and M. Painello, Phys. Rev. Lett 55, 2471 (1985).

[36] P. Ballone, W. Andreoni, R. Car, and M. Painello, Phys. Rev. Lett 60, 271

(1988).
[37] I. Stich, R. Car, and M. Painello, Phys. Rev. B 39, 4997 (1989).
[38] Pietro Ballone and Paolo Milani, Phys. Rev. B 42, 3201 (1990).

[39] A. Bulgac and D. Kusnezov, Phys. Rev. A 42, 5045 (1990); Phys. Lett. A
151, 122 (1990).

[40] D. Kusnezov, A. Bulgac and W. Bauer, Ann. of Phys. (204, 155 (1990).

 

[49

142
[4%

HM

M5]

[49

[47

L_a

14$

14%

501

[59

152

17
S. Bjornholm, J. Borggreen, O. Echt, K. Hansen, J. Pedersen, H. D. Ras—
mussen, Z. Phys. D 19, 47 (1991).

E. A. Rohlfing, D. M. Cox, and A. Caldor, J. Chem. Phys. 81, 3322 (1984).
R. E. Smalley, Laser chem. 2, 167 (1983).

H. W. Kroto, J. R. Heath, S. C. O’Brien, R. F. Curl, and R. E. Smalley.

Nature 318, 162 (1985).

T. P. Martin, T. Bergmann, H. Goehlich, and T. Lange, Chem. Phys. Lett.
172, 225 (1990).

H. Gohlich, T. Lange, T. Bergmann, T. P. Martin. Phys. Rev. Lett. 65, 748

(1990).

W. Kratschmer, L. D. Lamb, K. Fostiropoulos, and D. R. Huffman, Nature

347,354(1990)

R. J. Wilson, G. Meijer, D. S. Bethune, R. D. Johnson, D. D. Chambliss,
M. S. De Vries, H. E. Hunziker, and H. R. VVendt, Nature 348, 621 (1990).

J. L. Wragg, J. E. Chamberlain, H. W. White, W. Kratschmer, and D. R.

Huffman, Nature 348, 623 (1990).
S. Wang, and P. R. Busec, Chem. Phys. Lett. 182, 1 (1991).

P. A. Heiney, J. E. Fisher, A. R. McGhie, W. J. Romanow, A. M. Denenstein,
J. P. McCauley, Jr., A. B. Smith, III, and D. E. Cox, Phys. Rev. Lett. 66,

2911 (1991)

J.H. Weaver, J.L. Martins, T. Komeda, Y. Chen, T.R. Ohmo, G.H. Kroll, N.
Troullier, R.E. Haufler and RE. Smalley, Phys. Rev. Lett. 66, 1741 (1991).

 

18

[53] A. F. Hebard, M. J. Rosseinsky, R. C. Haddon, D. W. Murphy, S. H. Glarum.
T. T. M. Palstra, A. P. Ramirez, and A. R. Kortan, Nature 350. 600 (1991).

[54] M. J. Rosseinsky, A. P. Ramirez, S. H. Glarum, D. W. Murphy, R. C. Haci-
don, A. F. Hebard, T. T. M. Palstra, A. R. Kortan, S. M. Zahurak. and A.

V. Makhija, Phys. Rev. Lett. 66, 2830 (1991).

[55] K. Tanigaki, T. W. Eddesen, S. Saito, J. Mizuki, J. S. Tsai, Y. Kubo. and

S. Kuroshima, Nature 352, 222 (1991).

[56] Y. Wang, D. Tomanek, G. F. Bertsch, and R. S. Ruoff, Phys. Rev. B 47.
6711 (1993).

[57] G. Barton, and C. Eberlein. J. Chem. Phys. 95, 1512 (1991).
[58] G. Barton, Rep. Prog. Phys. 42, 963 (1978).
[5.9] A. Fetter, Ann. of Phys. 81, 367 (1973); ibid 88, 1 (1974).

[60] G. Bertsch, A. Bulgac, D. Tomanek, and Y. Wang, Phys. Rev. Lett. 67,

2690 (1991)].

[61] 1. V. Hertel, H. Steger, J. de Vries, B. Weiser, C. Menzel, B. Kamke and W.
Kamke, Phys. Rev. Lett. 68, 784 (1992).

[62] A. Herrmann, E. Schumacher, L. Wéste, J. Chem. Phys. 68, 2327 (1978).

[63] K. 1. Peterson, P. D. Dao, R. W. Farley, and A. W. Castleman, J. Chem.
Phys. 80, 1780 (1984).

[64] G. Delacretaz, J. D. Ganiere, R. Monot, and L. Wéste, J. Appl. Phys. B 29,

55 (1982).

 

19

[65] J. L. Gole, G. J. Green, S. A. Pace, and D. R. Preuss, J. Chem. Phys. 76,
2247 (1982).

[66] D. M. Lindsay, D. R. Herschbach, and A. L. Kwiram, Mol. Phys. 32, 1199
(1976).

[67] G. A. Thompson, F. Tischler, and D. M. Lindsay, J. Chem. Phys. 78, 5946
(1983).

[68] M. Hoffmann, S. Leutwuler, and W. Schulze, Chem. Phys. 40, 145 (1975).
[69] W. Ekardt, Phys. Rev. B 29, 1558 (1984); ibid 31, 636 (1985).
[70] M. Manninen, Phys. Rev. B 34, 6886 (1986).

[71] I. Moullet, J. L. Martins, F. Reuse and J. Buttet, Phys. Rev. Lett. 65, 476

(1990).
[72] U. Rothlisberger and W. Andreoni, J. Chem. Phys. 94, 8129 (1991).

[73] D. Tomanek, C. Sun, N. Sharma, and L. Wang, Phys. Rev. B 39, 5361
(1989).

[74] S. Sugano, Microcluster Physics, Springer—Verlag (1991).

[75] D. Tomanek, S. Mukherjee, and K. H. Bennemann, Phys. Rev. B 28, 665
(1983).

 

 

Chapter 2

Correlation Energy and Curvature
Energy of Nuclei

L 2. 1 Introduction

Direct adjustment of mass—formula to the ground state masses of nuclei through-
out the periodic table yields values of the curvature—energy coefficients ac (ac in the
term aCA1/3 in the mass—formula) close to zero. However many previous theoretical
calculations give a typical value 10 MeV for ac. Calculations within mean—field meth-
ods have yielded the value 9.34 MeV for the Seyler—Blanchard effective two—nucleon
interaction [1], values ranging from 9.52 to 12.99 MeV for six representative Skyme
interactions [2]. And also in a HF calculation of a sequence of large mirror nuclei with
the Coulomb force turned off, the value 11.86 MeV has been obtained [3]. Stocker et
a1 [4] have made a systematic study of the problem. From their analysis they con-
cluded that the mean—field approximations are unable to reproduce the entire set of
nuclear mass—formula coefficients when the saturation density and surface diffuseness
are constrained to their experimental values. However, there is one very important
effect, not yet considered, which might be responsible for this discrepancy; namely

the quantum corrections arising from the zero point oscillations of the nuclear surface.

20

 

 

 

21

It seems very natural to link the curvature corrections to the binding energy with the

surface oscillation, which modifies the local curvature when it is excited.

By summing up all the ring diagrams, we calculate the correlation energy Econ
for spherical nuclei in the RPA approximation with particular attention paid to the
A—dependence. We assume that the total ground state energy is the sum of Emean
and Econ, here Emean is the ground state energy in the mean—field approximation,
ECO" is the correlation energy in the ring RPA approximation. Since the fit of Em:
or other mean—field values of the total ground state energies gives the right values for
all the coefficients except aC which is far away from its experimental value zero, we
hope that the fit of E00,, to the same formula gives an c1C about the same magnitude
as that of the fit of Emean but with an opposite sign and all other coefficients to be
small compared with their corresponding values in the fit of Emean. In order to have
feWer coefficients to be fitted, we have used a set of imaginary nuclei with N = Z and
with Coulomb interaction turned off. Then we only have the volume term. surface
term, and curvature term to fit. The extracted curvature energy coefficient ac indeed

has the right order of magnitude and sign.

\

2.2 RPA Formalism of Correlation Energy

The formula to sum up all the ring diagrams is derived in many textbooks, e.g [5].
Because of the finite sizes of the nuclei, we will use the coordinate representation of

the RPA theory. In this representation, the formula for the correlation energy is given

by
1 ldA dw °°
Ecorr = _' __d3 d3 /_ AV On
22 o )1 :1: CE27rn:__:3( G)
1 1d/\ dw 1
: _- ___d3 3,,____/\ 03
220 A xdf2w(VG)1—AVGO

 

 

22

1
= 12 d—Adsa' (13.1: '21::

2 0 A (AVGO)/\VGRPA, (2.1)

with the RPA Green’s function given by

1

vRPA 0
G" :G 1— AVGO’

(2.2)

Note that the second order ring diagram is not included in the summation, since
it is already taken into account in the effective interaction in the usual mean—field
calculations. Also as it is usually done, we make an angular momenta decomposition
of GO and V, this procedure permits RPA to be calculated from independent radial
equations for each multipole, and also permits the ECO” to be calculated for each

multipole[7 ]. The following defines the multipolar interaction VL and Green 8 function

'0
GL,

V(r,r’) = ZVL(r,r’)YgM(f)l/},M(f’),
LM

000.71") = 2: G74“ TVS/EM“)YLM(7QI)/(Tr’)2'

Lil/I

Then the correlation energy is calculated as the sum of that of each multipole,

Lmar

Ecorr — _2( (2L +1)'EcLorr (23)

Here Lmax is the maximum multipole which is included in the summation, ECLON, is

the correlation energy corresponding to each multipole L and is given by

1d)1,dw

Egg, = 52 0 7.1m 2—(AVLC0 ) AVLG‘BPA, (2.4)
where
GIaPA_ GO 1 (2.5)

L1 — AVLGg

 

—7—_ -“-‘5-'- "'7 "1“"- ‘7 I "V

23

GEPA is usually found by matrix inversion if we treat both GL0“, r’) and VL(r, r') as
matrices in coordinate space. In the following, we will assume that the multipolar

interaction V1, is of separable form [6]

 

 

Vi(r.r')=a ‘ ' =af(r)f(r’). (2.6)

where v is the central part of the Wood—Saxon potential. The coefficient a is con-
strained by requiring GEPA to have a pole at w = 0 for L21, the spurious state. In
the case of separable interaction, 6?pr can be found analytically and therefore the
integration over A can be done analytically, this amounts to an enormous saving of

computational time.

Note that GRPA satisfies a Dyson—like equation

GR“ = G0 + GOVGRPA. (2.7)
For each multipole, we have

G’EPA(7~, r') = 016.7") + a / G%(vr,r1)f(r1)dr1 / d‘r2f('r2)G§PA(rg,r’). (2.8)

If we now multiply f(7‘) on both sides and integrate over r, we would get

. , _ f ali‘f(7‘)G0 (737")
/ ”(765109” ) “ 1—a ff(r1)G‘},(r1,Lr2)f(r2)dr1dr2' (2‘9)

 

Therefore 6'pr is given by

angU‘, 7‘1)f(7‘1)d7'1 fdr2f(r2)G%(T2,7‘l). (2.10)

RPA 7, 7,1 ____ 0 7' 7,!
GL, (’ ) GL( 7 )+ 1—aff(r1)G%(r1,r2)f(r2)dr1d7‘2

If isospin is included, GEPA, 0% and VL become 2 X 2 matrices. A similar argument

leads to

0’3“ -_.— G2+G},, (2.11)

24

 

G%(7‘,7") = ITO(7‘,7") ( (1) (1) ) , (2.12)
G'L(7‘,7") = 1_ Iguana/11% (7‘ 7‘1)f (7‘)1 )drI/drgf()()7‘21'l 0(7‘2,7")< ] ] ) , (2.13)
a = ff(7‘)II0(7‘,7")f('r’)d7‘cl7". (2.14)

Here H0(r, 7") is the free ph Green’s functions for protons and neutrons (in our case,
since N = Z and the Coulomb interaction is turned off, it is the same for both protons
and neutrons), so it can be factored out in front of the matrices. From the condition

that 6'pr has a pole at w = O for L = 1, we have
a =1/2a at w = 0 L =1. (2.15)

Finally the correlation energy is given by

1
EL" = T7‘ 44.6.6 —(/\VLG°)2 AVLGEPA
2i 0 /\ 271'
()Aaa3
42/()1(1 — 2/\—a—a 27r° (2'16)

Here T7‘ is the trace of the matrix. After the integration over /\ is done, the final

result is quite simple,
L 2 d“ 2
ECO“, = —z/-2——[2(aa) + 2(aa) + clog(1-— 2(aa)], (2.17)
77

clog is the complex natural logarithm. Of course the energy is real, and we should

take the real part of equation 2.17.

2.3 Numerical Details and Results

In our calculation, the single particle spectra and wave functions are obtained by

solving the Schrodinger equation using the Wood—Saxon potential. In order to be

25

able to treat the open—shell nuclei, we have included pairing and occupation numbers
for the open—shell levels in our calculation, the corresponding formulas for the Green’s
function are given by Migdal [8]. We have chosen the empirical value of the gap
parameter A = 1214‘”2 for all the openeshell levels, and A = 0 for the fully occupied
levels. The occupation numbers for the open~shell are calculated from the BCS

equation,

 

 

[\le—‘|

2_
Uk—

1— 0—6” . 1.. .
( \fléx-éflz‘l‘Ai) O 13)

The Fermi energy 6f is of course determined by

~2ng = N, (2.19)

k>0

N is the number of particles in the open shell and N = 0 for magic nuclei. The
separable potential has the same weight for all multipoles. But there must be a cutoff
at the high momenta. A natural cutoff for Lmax would be around kfR, this turns
out to be around 2141/3. Since 2A”3 is seldom an integer, we take the nearest integer
larger than 2A1/3. We further assumea Fermi—type weight function (1+ 6(L“L°)/L1)‘1
with L0 = 1.2141/3 for each multipole. For the integration over w, we have integrated
from w = 0 to w = 200 MeV with step size 0.5 MeV. The calculated energies and the
least‘square fitted coefficients along with their variance 02 of magic number nuclei
with different choice of L1 are shown in tables 2.1 and 2.2, and those of open—shell

nuclei are shown in tables 2.3 and 2.4.

2.4 Summary and Conclusion

The macroscopic part of standard nuclear mass formulas comes from an expansion

of nuclear energy on powers of (41/3. Such an approach is based on the geometry

 

26

of the nucleus, which consists of a bulk part separated from a thin skin, and is
not associated with approximations like HF, ETF, which however have been used
to reproduce empirical mass—formula coefficients. Direct adjustment of mass—formula
coefficients to ground state masses of nuclei yields values of the curvature coefficient ac
close to zero, whereas values of (1C resulting from mean—field calculations are always
close to 10 MeV. We have pointed out in this short chapter that Eco... has a very
different A—dependent behavior from that of the mean—field total ground state energy
Emean and that there is still room to account for the apparent discrepancy if higher

corrections to the mean—field calculation is taken into account.

 

Table 2.1: Correlation energies divided by 141/3 for magic number nuclei with

L1 = 1,2,3

Table 2.2: The least—square fitted coefficients of Econ. / 141/3 to the formula

27

 

E/A1/3(L1=1)

E/A1/3 (111:3)

 

 

ww
oooooll

50
82
126

 

-3.1120
-3.9128
-3.6554
—3.9785
-4.7833
6.0657

 

—2.9525
3.7988
-3.6105
-3.9963
-4.7949
-6.0786

 

-2.8785
3.7343
3.5724
43.9989
4.7881
-6.0864

 

ECON/Al/3 = aC + a3A1/3 + avA2/3, (1C is the curvature coefficient.

 

2

 

 

 

 

 

 

 

L1 aC as av a

1 -4.2535 0.8008 -0.1690 0.590E-1
2 -3.6783 0.5965 -0.1514 0.513E—1
3 -3.4686 0.5345 -0.1471 0.475E-1

 

 

Table 2.3: Correlation energies divided by 141/3 for open—shell nuclei with

28

 

 

L1 = 1, 2, 3
N=Z E/A1/3 (131:1) E/A1/3 (1.1:?) E/A1/3(L1=3)
18 8.1386 4.9301 —4.8263
22 8.1053 4.9685 4.8997
30 8.7112 8.4141 8.3198
48 8.1914 8.1599 8.1360
52 8.4939 8.4422 8.4045
80 8.8130 8.7859 8.7556
84 8.8491 8.8521 8.8565
122 8.7240 8.6967 8.6769
124 8.5299 8.5247 8.5203

 

 

 

 

 

Table 2.4: The least—square fitted coefficients of Eco." /A1/ 3 to the formula

EMT/Al” = 67C + asA1/3 + avA2/3, (1C is the curvature coefficient.

 

 

2

 

 

L1 ac as av a

1 -7.7925 1.4116 -0.l945 0.413E-1
2 -6.7973 1.0972 -0.1695 0.215E-l
3 -6.3635 0.9599 -0.l584 0.182E-1

 

 

 

 

 

 

Bibliography

[1] W. D. Myers and W. J. Swiatecki, Ann. of phys. 55, 395 (1969).

[2] M. Brack, C. Guet and H. Hakansson, Phys. Reports 123, 275 (1985).

[3] J. w. Negele and D. Vautherin, Nucl. Phys. A 207, 298 (1973).

[4] W. Stocker, J. Bartel, J. Nix and A. Sierk, Nucl. Phys. A 489, 252 (1988).

[5] A. L. Fetter and J. D. Walecka, Quantum theory of many—body systems, (San

Francisco, McGraw—hill, 1971)
[6] H. Esbensen and G. F. Bertsch, Phys. Rev. Lett. 25, 2257 (1984).
[7] G. F. Bertsch and S. F. Tsai, Phys. Reports 18, 125 (1975).

[8] A. B. Migdal, Theory of finite Fermi systems and application to atomic

nuclei, (New York, Interscience Publishers, 1967)

29

Chapter 3

Collective Electronic Excitations
in Fullerenes

The large number of active electrons in a carbon cluster and the strong Coulomb
interaction among them lead to a rich spectrum of collective states. States with total
angular momentum up to L 2 85 have a strong collective character. The equivalents
of both 77 and o plasmons in graphite are predicted for a whole range of carbon
clusters C20,60,70,100. The theoretical results compare very well with experimental

photoexcitation results in C60, 70 and EELS data on gas targets C6030.

3.1 Introduction

The discovery of C60 cluster and the successful synthesis of macroscopic quantities of
C60 crystal have spurred a flurry of activity in the study of properties of both C60
cluster and other fullerene carbon clusters [1]. Since these fullerenes are essentially
hollow graphite shells, they should support multipolar plasmon oscillations, closely
related to the plasmons in graphite. Indeed, the existence of the giant dipole plas-
mon mode in C60 was first predicted [2, 3] and then an independent photoionization

experiment on a C60 gas target [4] confirmed its presence.

30

 

31

The presence of the strong collective dipole plasmon in C60 suggests the exis-
tence of other multipolar excitations and the multipolar collective excitations in other
fullerenes as well. In this chapter we shall present an analysis of the collective exci—
tations in C20, C50, C70 and C100 and compare our predictions with available experi-
mental results on isolated fullerenes (gas targets). The dipole excitations are put into
evidence by photoionization experiments [4] and the the other multipole excitations

are evidenced by EELS experiments [5, 6, 7].

We should mention that there is a multitude of experimental results, concerning
plasmon excitation in pure solid C60 targets and doped solid C60 as well [8, 9, 10]. The
main difference between the plasmons in isolated fullerenes and fullerites [8, 9, 10]
and graphite and diamond [11, 12, 13] is mostly an upward shift of the o plasmon in
fullerites, by approximately 8 eV or more, because of the strong long range Coulomb
interaction between electrons belonging to different neighboring molecules. At the
same time, the properties of the 77 plasmon seem to be affected little by the presence

of other molecules.

3.2 Theoretical Description

We use the linear response theory (Random Phase Approximation — RPA), which is
the most appropriate theory for large systems with mobile electrons where screening
is significant, to describe the excitation of plasmon—like collective states in these
fullerenes. We determine the single-particle wave functions and energy levels using a
tight—binding (TB) Hamiltonian, which has been used to study the relative stability
of different carbon cluster structures [14] and the multipole response functions of

fullerenes [2, 5, 7].

 

32

3.2.1 Description of the Tight—binding Hamiltonian

The tight-binding Hamiltonian, which considers only the s and p valence electrons

of carbon, is given by [14]

H = anaiviaad + Z tag(r,-j)a:mag,j . (3.1)

0,! 07373.7].

Here, i labels the atomic sites and a = s,p;,.,py,pz labels the atomic orbitals. 60, is
the orbital energy, and tag are the hopping matrix elements between different sites.
The parameters have been obtained from a global fit to Local Density Approximation
(LDA) [15, 16] calculations of the electronic structure of C2, a graphite monolayer,
and bulk diamond, for different nearest—neighbor distances [14], similar to what had
been done previously to describe silicon clusters [17]. The diagonal elements of this
Hamiltonian are the energy levels 6, = —7.3 eV and cp = 0.0 eV. The off—diagonal
matrix elements tag(r) are assumed to have a distance dependence ~ r‘g. Their values
at r = 1.546 A, which is the equilibrium nearest—neighbor distance in diamond, are
V330 = —3.63 eV, Vspa = 4.20 eV, Vppa = 5.38 eV, and K,” = —2.24 eV in the
Slater—Koster parameterization [18]. In this Hamiltonian, we consider those atoms as
nearest neighbors which are closer than the cutoff distance rc = 1.67 A. This is the
average of the nearest— and second nearest—neighbor distances in bulk diamond, and

hence near the minimum of the radial distribution function.

The equilibrium nearest—neighbor distance between carbon atoms in C60 is 1.453
A on pentagons (“single bonds”) and 1.369 A between pentagons (“double bonds”)
[2], corresponding to a buckyball of radius R z 3.5 A. The spectrum of C60 obtained
using this tight—binding Hamiltonian is shown in fig. 3.1. This spectrum is in good
agreement with previous ab initio calculations. The total width of the occupied band
is 19.1 eV, to be compared with the LDA values of 18.8 eV of ref. [19]. The level

ordering near the Fermi level agrees with results based on the LDA and other methods

 

2325

m

V
5.

tlr

(

10

(LUMO)
(HOMO)

 

E(eV)

 

 

 

3 ”ll" II II llllll

‘0 II ll

0
I
‘Q

2’

J‘
J’

T!

'3 (LUMO)

 

   

o 00-00-000000..'_— ”u‘HOMO)

.c...mooo_ n

-37- fl .

 

Figure 3.1: (a) Tight—binding Single-particle energy level spectrum of a C60 cluster.
The levels have been sorted by symmetry. (b) Expanded region of the energy level
spectrum near the Fermi level. Allowed dipole transitions hetween states with gerade
(g) and ungerade (u) parity are shown by arrows [From C. FBertsch, A. Bulgac,
D. Tomanek, and Y. Wang, Phys. Rev. Lett. 67, 2690 (1991)].

 

34

[19, 20, 21, 8, 23]. The gap 2.2 eV between the highest occupied (HOMO) and lowest
unoccupied (LOMO) level compares well with the LDA values of 1.8 eV from [19] and

the experimental value [8] of 1.9 eV.

The HOMO to LOMO transition is forbidden by parity, and the lowest optically
allowed transitions are Hu —> T19, Hg —+ Tm, and H, —-> Hg, with the tight—binding
excitations energies of 2.8 eV, 3.1 eV, and 4.3 eV. These results compare well with

the LDA values 2.9 eV, 3.1 eV, and 4.1 eV [19].

3.2.2 Linear Response Theory of the Tight—binding Hamil-
tonian

The free particle—hole propagator in the energy representation is defined by

Z < rlph > 2(ep — 51,) < phlr' >
p. (a. —54)2 — (w + 27272

 

. (3-2)

Go(r,r’;w) =

where p and h denote particle and hole states, spy, their corresponding energies, to

the excitation energy and 7] is an imaginary part.

The TB single—particle wave functions have the following structure

1
$7741.): 22‘: fn(R i)\/4_7T¢O 0(T1) +gn(R M\/:¢1() i)] 7 (33)

where n denotes either a hole state or a particle state, 1‘ is the electron coordinate,
R, are the carbon ion positions, r,- = r — R,, the summation is over the carbon sites,
050,1(7‘) are the LDA carbon 2s—~ and 2p— electron wave functions. fn(R,) and gn(R,-)
are the amplitudes to find an electron at site i either in the 2s—state or in one of the

three 2p—states. The particle—hole state then reads

< PIPh >= fi: [fp(R4)fh(R.-)¢3(r7) + 87434) ° gh(R4)¢f(r.-)] YooU‘i) +

47r Lif— z‘).gh —m( Rt) + fh(Ri)gp,—m(Ri)l ¢o(7‘r)¢1(7‘i)Ylm(f‘r) +

 

35

3 A
107 ZHY” [gm-7 >< 84111012,-.. 81471484.). (3.47
In the above formulas we have used interchangeably Cartesian and Spheri cal notations.

The Random Phase Approximation (RPA) Green’s function is determined as the
solution of the integral equation G = Go—GOVG, where V = e2/|r—r'] is the Coulomb
interaction among electrons. The matrix elements of the Coulomb interaction are

computed through its Fourier representation

2

e (‘32 1 .
= [r4 — 1‘: + Rk -- R1] : 27r7/d3q32—8XP12CI' (rk — 17+ Rk - RH], (3'5)

 

V

where rkJ are the electron coordinates with respect to sites 10,1 and R1,; are the
coordinates of the corresponding carbon ions. The response function (transition
strength) of the system to a weak external one—particle multipolar field F(r) is given

by S = Im < FlGlF > /7r.

In computing the matrix elements of the residual Coulomb interaction arising from
the l = 0 transitions in Ref. [2] for the dipole response, it was assumed that the spatial
extent of the s— and p—states is vanishingly small, this approximation results in a
slight overestimate of the Coulomb matrix elements. Another approximation in ref.
[2] was to neglect the l = 2 component in the particle—hole state [ph > (see eq. (3.4))
resulting from a p—hole state to a p—particle state transition. These approximations
make their approach unsuitable to study other multipolar responses and also the
multipolar responses of other nonspherical fullerenes. The present calculations are

performed without these approximations.

Since only dipole excitation is optically allowed, one can put into evidence higher
multipole plasmons in an inelastic electron scattering experiment. In the Born ap-

proximation the differential cross section of electron scattering is [24]

d2” 62m 2 4P, - 2
dfldw 2 7,2“ 7532' <n|exp(—zq-r)|0> | awn—130w)

 

 

36
82m 2 4p’
2 -— — a), , 3.6
where p and p, are the initial and final linear momenta of the electron, q = p — p,
is the transferred linear momentum, m the mass of the electron, w the energy of the
excited state and S(w,q) the spectral function of the cluster, which depends solely

on the properties of the cluster.

In the following two subsections, we will outline the application of the RPA theory
of a tight—binding Hamiltonian to a spherical cluster (C60) and nonspherical clusters
(C20, 70, 100). The detailed treatment of C60 is given in Appendix A, and that of other

nonspherical clusters in Appendix B.

3.2.3 Outlines of the Calculations of C60

The treatment of C60 is relatively simple. Because C60 has nearly perfect spherical
symmetry, good angular momentum quantum numbers can be introduced [2, 5, 7],
and the the RPA equation is solved for each multipole (see Appendix A). The total
angular momentum of a ph—state (see eq. (3.4)) has two components L and l, i.e.
J = L + l, where L comes from the fn(R,-) and gn(R,-) amplitudes alone and 1 from
the transitions on each site: I = 0 from s —> s and p —> p; l = 1 from s —> p and p —> s;
l = 2 from p —-> p. The three terms of the pit—wave function correspond to l = 0, l = 1
and l = 2 transitions on each site. Since we will be concerned with natural parity
states only and the spin variables merely double the number of available states, the
following selection rule applies (-1)J+L+l = 1. As a result for J 2 2 only the following
(L,l) configurations are allowed: (J, 0)2,(J — 1,1),(J + 1,1),(J — 2,2),(J, 2) and
(J+2,2). For J = 0: (0,0)2, (1, 1) and (2,2) and for J =1: (1,0)2, (0, 1), (2, 1), (1,2)
and (3,2) are allowed. The two (J,0) configurations arise from two I = 0 types of

transitions, see eq. (3.4).

37
Using the expansion of a plane wave into spherical harmonics, we can rewrite the

Coulomb interaction eq. (3.5) as

2 2

 

e e A
V— In —r2+R1 _R2|_ —2 2/exp[ik(r1—r2+R1 —R2)]dkdk
(471')3e2
: _27r_2—_/LzlL jL (kR1)YL]W( (R1)YL1\/[(k )Z Z- LjLI( ()kR2 YLI\JI(R2)
L’M’
YL’M’O‘ )ZZJ1( ((kT1)Yzm (1‘1) MEI-l, Jz'f (INDY I’m ,f‘( 2)}i'm'(f<)dkdf9 (3-7)
[m l’m’

After the integration over k, PM is done by using the ph-wave function eq. (3.4). one
can represent the residual Coulomb interaction as a matrix in the representation (L, I)
described above. In a similar fashion to Ref. [2], the integral equation G = Go—Go VG
is then reduced to a matrix equation in the (L, l) representation of dimension 4 x 4
for J = 0, 6 X 6 for J = 1, and 7 X 7 for J 2 2. The resulting matrix equation
is solved by matrix inversion The multipolar responses are obtained by applying the

multipolar fields F(r) = r’ng(f‘) (F(r) = r2Y00(f‘) for l = 0).

In the case of spectral function, the external field is the plane wave exp(iq - r),
here q = p — p, is momenta transfer. It is computed by expanding the plane wave
into spherical harmonics and all excited states with angular momentum up to L = 20
are included (see Appendix A for details). We have checked that the responses and
the spectral function of C60 computed either by assuming spherical symmetry (as
described here and the details in Appendix A) or by taking into account the actual
3—dimensional ionic configuration of the C60 molecule (as described below and the

details in the Appendix B), lead to essentially the same results.

3.2.4 Outlines of the Calculations of C20,70,100

Since these clusters do not have spherical symmetry, good angular momentum quan-

tum numbers can not be introduced for them, and because of that the actual numerical

38

implementation of the RPA is computationally more demanding (e.g. the calculation
of the whole dynamical structure function for C60 requires 1—2 hours, while performed
in the 3—dimensional manner described here and Appendix B, it requires about 300
hours on a CONVEX—240). To this end, we have formulated the RPA method for ar-
bitrary 3—dimensional structures in the TB approximation and used it to compute (or
recompute in the case of C60) the electromagnetic response of isolated fullerenes (the
details are given in Appendix B). The matrix elements of the Coulomb interaction
are computed through its Fourier representation, but this time instead of treating R1
and R2 as two different vectors (as we do in the case of C60), we treat R1 —- R; as

one vector, and perform the expansion of the plane wave into spherical harmonics

 

 

accordingly,
(,2 62 .
v = In __ m + R1 __ R2] = fi/expfik - (r1 — r2 + R.1 — R2)]dlcdk
(473382 .,. . , ~ ._,,. , . .
= 2.2 / gzyz<knmm(rlmm(k)§z JI'(k7"2)Yz'mI(I‘2)Yl'm'(k)
LEA/:1iLijkR12)YLM(R12)YEM(R)dkdf( (3-8)

Instead of solving the RPA equation for each multipole, we solve the whole RPA
Green’s function C(r, I") as a function of r and r’. Essentially, the integral equation
is transformed to a matrix equation in the (lmi) basis after the integration over k, rim;
is done , here I = 0, 1, 2 denotes the particle—hole transitions on each site i and m
denotes the corresponding third components of I. There are two terms for l = 0 (see
eq. (3.4)), three terms for l = 1 and five terms for l = 2. Therefore there are 10
lm terms for each site i and 10N terms in the particle-hole states in the (lmi) basis,

where N is the number of atoms in the cluster. A particular multipolar response 1 is

 

1strictly speaking there is no one—to—one correspondence between the multipolarity of the external
field and the corresponding one for the excited plasmons, the results presented were averaged over
all possible spatial orientations.

39

obtained by applying a multipolar external field. The spectral function is computed
by applying the external field exp(iq- r) and averaged over the direction of q, see

Appendix B for details.

3.3 Theoretical Results and Comparison with Ex-
perimental Data

In figs. 3.2—3.5, we present the computed free (no residual Coulomb interaction among
electrons) and the RPA responses of C20, 60, 70, 100 clusters to a multipole external field
F(r) = r’Y1m(f‘), for I = 0 —— 8 (F(r) = 7‘2 for l = 0). C20 has only 12 pentagons,
while C70 and C100 can be obtained basically from a C60 molecule, by cutting it into
two equal semispheres and joining them by one and four additional carbon rings
respectively. We compare our predictions of the photoexcitation probabilities of C60
and C70 with the measured photoionization spectra [4] in figs. 3.6 and 3.7. One
should note the good agreement between our predictions and the experimental data,

even the width seems to be well obtained.

It is obvious from the these figures that a strong screening is obtained at low
energies for all of the multipole responses of all these carbon clusters. Except for the
monopole responses in C60 and C70, all other multipole responses have a very similar
structure, a low—energy mode around 6-10 eV and a high—energy mode around 18—22
eV. These two modes are the obvious analogues of the 7r and o plasmons in graphite
[11, 12, 13] and the corresponding plasmons seen in solid C60 [8, 9, 10]. In both
graphite and solid C60 the o plasmon is at a higher energy than it is in an isolated

fullerene.

The states with higher angular momenta show no collective behavior or in other

words the corresponding multipole responses are single-particle in character. The

 

~10

multipole responses of C30

 

 

 

 

 

 

 
 

 

 

 

 
 
 

  

      

 

 

 

 

A
“ """"""'f"vvfiff'frvvr"" 'VVTI""V'V""'
”g ’1 M 1 L-1 1 m
' i 7‘ i i.
'9 - \ I
3 ’ i '1
7 I | l
g .' , - \‘ u .. 1,, r-
fl 1-3 11.4 [ L-s [
t 7.1 1 1.
III ’ | I
3 [.lflf‘ 1“ 1|! 1‘
i i. “ 4 'I ""
a L 1 m | W 1 m
3 7‘ 1 II
8 _:\l\. 0 a 1‘
g {\IAAIAAAALA-AA --nm---1-nnn LEA----ILA-J -
O 10 20 30 0 10 20 30 0 10 0 30

excitation energy (eV)

Figure 3.2: Multipolar free (dashed lines) and RPA (solid lines) responses of C20

41

multipole responses of Coo

 

   

 

 

. "V'fi7"'f'f'U"""[""‘f"'I'
L 1 Lil I LIZ
D
D
t
: 11
I
: 1| "1
D
._ l\ 741 \
E f r- I / \
7 \
i' “H \ DIG

 

 

 

 

   
 

 

 

 

 

  

1AA4-1; Am 44 I; A l

o “10 20‘ 300 210200800 ‘10 2003230
excitation energy (eV)

diflerential oscillator strength (arb. unit)

 

Figure 3.3: Multipolar free. (dashed lines) and RPA (solid lines) responses of C360

 

multipole responses of C70

 

 

 

 

 

 

 

 

 

 

 

 

9

.3 r r tho I 1:.1 t I vh‘vvrvvi-vyvvv'v

. " 7‘

1° I \ ’1

h I 1 1

v I \ 'v I ~

g r, i\\ H ’ ‘ M
,1 L-S ,\L-4 ,\ M

‘3 7‘ 7‘ 7‘

u 7‘ 1‘ ‘

a I \\ I \ r \

8 ,\ L-e ,1L-7 ,\ m

5 7 ‘ 7 7
7 _ 7 I

g ’--.----.L\.i ’12“--. 1: ’U- ---.: -
0 10 20' 30 O 10 20 30 O 10 20 30

excitation energy (eV)

Figure 3.4: Multipolar free (dashed lines) and RPA (solid lines) responses of C70

 

multipole responses of C100

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

09 fl r""rffivri*'* ' W"*I*" r' vvvvvvvvvvvvv
a :1 M 7 LI! 1 we
‘ 7‘ 7.,
' :1 \
a I." ‘s i 1 I‘
v i." / i i \
g E. / \\._ / \\ l JR \
. - 7‘ us 7‘1.“ 71 w
a " H 11
h "1 ’\ l \
3 .1 ,- 1. \ p I \«‘._~
8 7‘ M [\L-7 I\ no
3 7\ 7\ 7\
7 \ 7 \ \
5 I -- 7 -4 7
g ’xxl---xlx§-~x . Ixml+-n+:im\mx [xxlxxAxlx\x>A A
O 10 20 300 10 20 300 10 20 30

excitation energy (eV)

Figure 3.5: Multipolar free (dashed lines) and RPA (solid lines) responses of C100

(orb. unit)

~44

 

  

 

 

err" - . Wes
l

D

> 0"“.

5 l,- . 75...... .' ' C60
?

C

4r

P

P

D

3r-

D

b

l

2L

p

r

l

1:-

b
ob41-44-1¥--Al----LAJALlA

10 15 20 25 30

excitation energy (eV)

Figure 3.6: Dipole RPA photoexcitation probability (solid line) and measured pho—
toionization cross section (arbitrary units) [4] of (7360

(orb. unit)

 

 

10 15 20 25 30
excitation energy (eV)

Figure 3.7: Dipole RPA photoexcitation probability (solid line) and measured pho-
toionization cross section (arbitrary units) [4] of (370

 

46

physical reason is relatively simple. For these multipoles the wave length of the
excitation is comparable with the C—C bond length and no electron collectivity exists
at such length scales. Take, for example, C60, the maximum expected multipole of a
plasmon Lmax = qmarR = R7r/d % 8, where R is the radius of the buckyball, d the
C—C bond length and qmax = 7r/d is the momentum of the excitation with a half—wave
length equal to the C—C bond, compares extremely well with the RPA results in fig.

3.3.

As remarked earlier, one should note that strictly speaking only in C60 there is
an essentially one—to—one correspondence between the multipolarity of the external
field and the corresponding one for the excited plasmons. For other fullerenes, which
lack spherical symmetry, the presented results were averaged over all possible spatial
orientations of the cluster with respect to the laboratory reference frame (or in other

words the orientation of the external field).

Since C100 is a rather deformed object and one might expect that one can detect
the deformation by observing the splitting of the dipole plasmon. However, as our
results show, when one applies the external field either along or perpendicular to the
symmetry axis of this cluster, the dipole response changes insignificantly, see figs. 3.8,
3.9, in contradistinction to what one sees for example in deformed nuclei or metallic
clusters [25, '26]. To some extent this is surprising. This feature of the plasmon modes
in fullerenes simply reflects the specific character of the structure of the single—particle
wave functions, which are strongly localized around the ionic shell. The curvature
of the shell is less important than the 7:" or a character of the single particle wave

functions, and therefore of the corresponding plasmon modes.

Recently [27], the energies and oscillator strengths for a number of allowed dipole

transitions have been measured for C60 in solution. A comparison between the calcu-

9
on

P
to

.0
H

O

differential oscillator strength (arb. unit)
. O
'0

Figure 3.8: RPA (solid line) and free (claslu-xl lino) dipole response of (‘100 when tllt‘

.°
a

P
.p

to

dipole

47

response along the symmetry axis

 

 

 

l

i 0100
D l‘

: \

: I\

. I \ ‘

D I \

T i“; \ \

I l \ //

f l \ / \

. l \

I:

r! a

,1 ~-
;/

l L L n 4 A n L L L
0 10 20 30

excitation energy (eV)

external field is along the symmetry axis.

 

dipole response perpendicular to the symmetry axis

 

 

 

 

7.?

a 006 ’ f v ' ' t ' ' __ v y - v v t v v

g 005 . Cl” .:
i / ‘ (

V I \ (
r / \ .

5 o 4 ’. I :

g . i I \ I, \\ 4

0 I I \

5 r I i ’ \ .

°° 0.3: ,\ I \ I ‘ \ ‘

h , \ I \ l

3 ’ l l \ I \ *

g . I \ \ l

a 0.2 f g \ ll \ .

8 : ' ’ \x

O , l \ \

s » I ‘ ‘ ~

5 , ‘

h 0.0 g l A A_ L; A L L 4 x A l g

.2 O 10 20 30

'3 excitation energy (eV)

Figure 3.9: RPA (solid line) and free (dashed line) dipole response of (7.00 when the
external field is perpendicular to the symmetry axis.

49

Table 3.1: The free and RPA dipole response of C60 (energies and oscillator strengths),
computed using two different tight—binding (TB) parameterizations [14, 29] compared
to the calculations of Ref. [28] and experimental results [27]. A * next to a transition
means that several close in energy levels have been represented as one. A + sign
means that the oscillator strengths of the two or three consecutive levels have been
added. The measured oscillator strength for the 3.04 eV transition is 0.015i0.005.
The last row is the sum of the oscillator strengths quoted in the table.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Free TB [14] [[ RPA TB [14] [[ Free TB [29] [[ RPA TB [29] [[ Ref. [28] Ref. [27]
[[eV [7 [ev [f [[ev [f [[ev if [[ev |f [[ev [f
2.88 4.04 2.94 0.027 2.6 3.5 2.68 0.015 3.4 0.08 3.04 0.015’
3.12 6.09 3.74 0.66 2.82 5.41 3.37 0.43 4.06 0.41 3.30
4.28 3.95 4.78 1.37 3.84 4.01 4.39 1.24 4.38 2.37 3.78 037
520+ 522+ 4.83 0.029 470+ 406+
526+ 528+ 0.64 4.93 2.70 5.17 0.44 507+ 0.30 435+ 010
537+ 2.88 5.69 1.15 538* 1.36 5.32 0.13 5.24 7.88 4.84 2.27
5.86 5.94 0.21 6.03 0.35 5.47 0.12 5.54 0.22 5.46 0.22
597+ 0.86 6.13 2.48 6.31 0.24 5.80* 2.86 5.78 10.74 5.88 3.09
6.73 6.26 0.50 6.14 1.64 628 6.36
6.76 0.40 6.40 1.99
6.90 0.15
|[ [ 17.82 [[ [ 7.49 [| [ 17.57 [[ [ 8.90 [[ J 22.00 [| [ 6.07:

 

 

50

lated free and RPA predictions, using two different parameterizations for the tight—
binding Hamiltonian, the calculations of Ref. [28] and experimental results of Ref.
[27] are presented in the Table 3.1. Even though a detailed correspondence between
computed and measured transitions cannot be established yet, the significant dy-
namical screening emerging from the RPA description compares very well with the
observed transition strengths, including the large suppression of the first allowed
dipole transition. It is unclear to us however, to what extent the computed prop-
erties of an isolated C60 cluster are related to absorption spectra of C60 in solution.
A measurement of the gas form will be highly desirable. The static polarizability of
a byckyball, OFT-63 = 195 and 207 and (13122) = 56 and 57 using parameterizations

[14, 29] respectively, shows the same order of rather strong screening.

The dynamical structure factors for C60 and C70 are represented in figs. 3.10 and
3.11. One can see that they are very similar, there are two “mountains”, one at an
excitation energy around 6—10 eV and a second one around 18—22 eV, corresponding
to the 7r and a plasmons respectively and their corresponding energies increase mono-
tonically with the increase of the transferred linear momentum. The profile of the
spectral function in the q direction is rather structureless, due to the fact that states
with different angular momenta are very close in energy. The dispersion laws for the
71' and a plasmons, inferred from the profile of 3(w, q) (u),r z 6 + q2 and 020 z 18 + 2q,
when energy and momentum is in eV and A‘1 respectively), are rather distinct
from classical estimates for a charged spherical shell [31]. The maxima of the spectral
function are at transferred momenta z 1.5—2 A“, indicative at the fact that states
with L % qR R: 5—7 are mainly excited in this process. The a plasmon always has
a higher intensity. At relatively high transferred linear momenta we do not expect

these computed dynamical structure factors to be very reliable though, because of

51

the limited TB Hilbert space we have used, but on the other hand we do not expect
a significant collectivization of the transition strength either when the transferred

momentum becomes comparable with Tr/d, where d is the C—C bond length.

In figs. 3.12 and 3.13, we compare the RPA predictions for the excitation of the
plasmon modes in C60 and C70 with recent EELS experimental results on a gas target
at an incident electron energy of 1 keV and at 1.5 deg for C60 and 1 keV and at 2.0

deg for C70 [5, 6, 7] .

The linear transferred momentum in the C60 experiment is "~V 0.41 A”, this corre-
sponds to an almost equal contribution from the dipole and quadrupole excitations,
and the linear transferred momentum in the C70 experiment is z 0.55 A”, and it cor-
responds mainly to the quadrupole excitation. But in both experiments, the excitated

multipoles are far from the main maxima of the spectral functions.

A few words about the experimental setup are in order. The detector used in both
experiments has a circular finite area, it forms a cone with the scattering center of
the sample, the surfaces of the cones are 1° off their axis. The axis of the cone in
the C60 experiment is l.5° off the beam of the incident electrons, and that in the C70
experiment is 2.0° off the beam of the incident electrons. The actual geometry of the

detector has been accounted for in our calculations for the differential cross sections.

In figs. 3.12 and 3.13 the computed differential cross sections of the plasmon ex-
citation are compared with the EELS data on C60 and C70. Since only relative cross
sections were available and the experimental data have been thus renormalized. One
can see that both the shape and the position of the plasmon modes agree very well
with the experimental data. Moreover, even some rather fine details of the experi—
mental results seem to be reproduced by the theoretical model. In photoexcitation

experiments one can excite only dipole states, while with electrons one can in principle

 

 
  

 

  
 
     
    
           
    

    

  

4 v v v v I f f r v r v v Via I v v v ' v v v v ' v
p O 0
g . 0 ..0
0 g . g I
l I I : I
O 0 . .00000'
t ..l I. .I' "I.. I
I 'I IIIII 'I
b .0 g. I... ...I. .0
I. I. .0 .. ....oo.... ..
3. .. .. . 0' .0......0.... '0'...
I. 'I . ' 0'. .'I
I .0 g. g
3 I
I I
I I a I
I I I.
I I
I. I I. I
A f . 0. .. 0
on I' 'I I.
p O o 'I
| .. . .0
I II
‘ ' I ' .'I
II . I. .. I
v . 0 0. .II
’ I ° 'I I"'
I g '0.
O‘ ' I ' '0.
I .00
I
' I ' "I
D . O. . I
. 0 g. .0.
i . . 0.
I
. I I. I
1 P 0 :.. . I...
. 'I III...
, . 0 0. ..'..0 0
' . " " I. III I
l O. I. ... II .... I. . 0.... ....I.. ..
I. .l I ..III.. I. III.. II I' I II
I .0 .0...l:...l... ..'.0.............I.I.I.I.O....'...I.......... 00......
.. 0...0....0000.0.... .000000.....'.'00.. ....0..
P 'IIIII" "IIIIIIIIIIIIIIIIIIIIO'

 

o A A L A 1 A A A A L A A A A 1 A A A A 1 A A A A L A

 

 

 

0 5 10 15 20 25
0: (eV)

Figure 3.10: The dynamical structure factor (spectral function) of (360

 

the dynamical structure factor of C"

 

   

 

 

4 - - - + . - . - - r -.
' I .I ...IIIIII....
. ° ' I' .I'
I
’ IIIIII.
.. '. 0.....|.OIO'OOII..
f .0 . .0 .............0
b I ' ' ....IIII|..
. I
3. I '
. I
I .0 ...
l I ' °
I ' '
A i I ‘ 'I
I ' 'I
CI. . I I.
| i I ' 'I
I ' 'I
‘ z. I ' 'I
V . g 0'
i I ' 'I
0‘ b ' ' 'I
I ' .0
, I ' :I
I
i o 'i
f I ‘ :I
1' .' -' '1
I ' 'I
i I ' 'I
0 . .gI
' I
I. ....I.I
i ..OIII...'.'.O
D
A l

 

 

, Figure 3.l l: The dynamical structure factor (spectral function) of ('70

54

excite any multipolarity, depending on the amount of linear momentum transferred
to the target. In the semiclassical approximation, the angular momentum of the ex-
cited state is given by L z qR, Where L is the multipolarity of excited state, R is
the radius of the target and q = pf,n — pin is the transferred linear momentum. The
EELS C60 data in fig. 3.12 thus correspond to approximately an equal contribution
of dipole and quadrupole states and the EELS C70 data in fig. 3.13 correspond to
approximately quadrupole state, therefore these are the first experimental evidence

of quadrupole states in isolated C60 and C70 molecules.

3.4 Discussions

In an RPA treatment one cannot describe the fragmentation of the oscillator strength
due to the coupling to more complicated states. The Landau damping (i.e. fragmen-
tation of the strength due to the 1 particle - 1 hole background) is accounted for
in RPA, but this mechanism is not sufficient to explain the significant width espe-
cially of the 0 plasmon, see the results presented in Ref. [2]. Since these clusters
are essentially 2—dimensional objects, the damping due to more complicated particle
— hole configurations was imitated with an imaginary part 77 z 02/8 eV in eq. (3.2)
(chosen as to describe the experimentally observed width of the dipole plasmon in
C60), characteristic of the coupling of the RPA modes to surface electronic oscillations

only [32].

However, one cannot fail to observe that the predicted peaks are about 2—3 eV
lower than the measured ones. We suspect that the reason for this discrepancy
lies in the very special structure of the plasmon spectra of these clusters, namely
that depending on the amount of linear momentum transferred to the cluster by the

external field, 202.,r and 37.07r are very close in energy to 020, see the dynamical spectral

 

 

 

-.,fl....,.-..r...l..fi,.
mo; , -
: 3:. C” ‘
L .9 '. ‘
80b .‘I. -
. i
B. " i I i
> 60 - " ' i
o - ' .. ‘
..\ C .' ‘1 -. «:4 o'.~ 3“.“ - 3
3 * i .0 M" ' 2"".1' ' '3'!“ .. ‘
c: ’ " .a-cé"'~".' "'- ' '- -'o' ‘
@40- ' .'.f' ' "fl -
3 ’ .s 3, 4.? 1?. *
b ..W 4
\ I. o I" 0' .

“b t figs.“ :
vzoe V"
,. . I.

Figure 3.12: Computed (solid line, absolute values) and measured (points) inelastic
electron excitation cross section of plasmon in (-350 with incident 1 keV electrons at
1.5 :1: 1°. In the computed results, the actual geometry of the detector has been
accounted for.

 

7) (i

 

 

 

dad/dado (Ag/eV.“)

60 P t v 1' I r i r I f I v I T v v I’ I l r v v r I r— ‘
H.
I 3%”: 070 I
- I ’ .—
50 P I Q 00' . O . . ‘
b ’ . ' . . . . . . 1
. 0° .' . ' ‘ "i, o o ‘
D ’ . I ‘fifi'z. . aitv.. ‘0. .‘ ‘
4O "' . .0. . ”5‘52: " . ' . 9 gr: —4
. V 'u’:’.' o "'0 . . ' ‘Pfi' ‘
p '. o o h 0‘: h‘ 0" .. '0. 00 1
P 0. o g ":0 (0".0‘3... . ‘ “‘0' ‘
30 '_ : .2 @319? ': “i" J
. ; '7; .3 ' 1.4;. i
b g. ' . '0 1
L I ‘92,? Wt:-
20 _ ' '1": '«t
I i?
b ’ .0
10 f 4
. 1 . .i

 

 

O A L l A A A A L A 1 L A L L L L L [A L 1

5 10 15 20 25
a) (eV)

Figure 3.13: Computed (solid line, absolute values) and measured (points) inelastic
electron excitation cross section of plasmon in ('70 with incident 1 keV electrons at
2.0 i 1°. In the computed results, the actual geometry of the detector has been

accounted for.

 

57

functions in figs. 3.10 and 3.11. In this respect it seems that these clusters are
very peculiar and the structure of the plasmon spectra indicates that there may be a
strong, essentially resonant coupling between these modes (namely 2w,r or 3w, with
too), unlike other many-fermion systems studied until now. A plasmon is similar to
a photon in QED. In QED it is known that the Furry’s theorem forbids the 7 ——> ‘27
process (7 is a photon), but at the same time the ’7 ——> 3’7 process is allowed. For
similar reasons (in particular because there is an approximate symmetry between
the particle and hole spectrum) we expect that a —> 2 7r process should be to some
extent less important than the a —> 3 7r process, see fig. 3.14. The finite geometry
of the carbon clusters, phase space considerations and the fact that a plasmon is
not exactly a free photon may lead to the fact that these two processes are actually
equally important. The coupling between the lp—lh and 2p—2h and 3p—3h electronic
excitations may also lead to a rather strong non—linear optical response of the carbon
clusters and to an upward shift in energy of the high—lying mode, besides the large
fragmentation of the strength of the a—plasmon. A proper theoretical treatment of
these effects goes well beyond the framework of the RPA and it is a rather involved

computational procedure (never before attempted in the literature to our knowledge).

A couple of aspects of our theoretical treatment deserve a few comments. In the
present approach we have approximated the continuum states with a set of discrete
states and thus we are not able to describe electron escape widths. On the other hand
one can suspect that a proper treatment of the continuum might change significantly
our results. We do not expect this to be the case. In Ref. [2] the continuum has been
accounted for in the jellium approximation and the electron escape widths turned out
to be much smaller than the fragmentation width of the plasmon. We have also ne-

glected exchange—correlation contributions to the residual Coulomb interaction, which

 

 

 

 

58

are much smaller than the direct Coulomb part, which has been taken exactly into
account in the present approach. Yabana and Bertsch [33] have analysed explicitly
the role of certain Coulomb exchange diagrams in inelastic electron excitation of C60
and found that they are relatively unimportant, even at small energies of the incident
electrons, where one might expect the exchange Feynman diagrams to be relatively

important.

The most debatable approximation is the use of a limited TB space for the elec-
trons. As a consequence the total oscillator strength accounted for is only a fraction
of the total (one has about 70 instead of 240 for C60 and approximately the same
ratio is obtained for the other fullerenes). This is a rather serious drawback of the
present approach. However, we have reasons to expect that the missing strength
is at much higher energies, than the ones considered here. First of all, the present
formalism seems to describe quantitatively pretty well the oscillator strength in the
low energy region, where absolute measurements are available [2, 5, 27]. At the same
time the relative intensity of the 71' and a plasmons (as seen in EELS experiments on
gas targets, see Refs. [5, 6, 7]. and the present results, where only relative intensities
are available so far) seems to be reproduced also very well. In a recent study of the
plasmon states in C60 in the framework of the jellium model, including however the
icosahedral distortion of the single particle spectrum, where most of the oscillator
strength is accounted for in the formalism, Yabana and Bertsch [34] have shown that
in the RPA up to an excitation energy of approximately 30 eV only about 30 ‘70 of the
total energy weighted sum rule, i.e. approximately 80 units of oscillator strength out
of 240, are present. The rest are located at much higher excitation energies. There
are some experimental indications as well that this might be the case, from the double

(and even triple) photoionization experiments on (360/70 gas targets [35]. A theoreti-

 

5 5)

11' 11'
0’ 0' 1T
1? 1I'

Figure 3.14: Feynman many—body diagrams responsible for the fragmentation width
of the plasmon modes in carbon clusters. Solid lines represent the electron propaga-
tors and wavy lines the plasmon modes, while the residual interaction is represented
through a point vertex.

(it)

Figure 3.15: Feynman diagrams of the decay of a plasmon through the emission of 2
or 3 electrons

61

cal analysis of the double and triple photoionization processes is warranted. One can
expect that a highly excited plasmon can easily decay through the emission of '2 or
even 3 electrons, in a 2— or 3—step process, see fig. 3.15. When a high energy plasmon
(more than 30 eV) decays through the emission of an electron (particle state), the
residue is left in a relatively deep hole state or again in a relatively highly excited
particle state, which subsequently might easily couple to a lower energy plasmon,
which subsequently decays through the emission of a second electron. This type of
process seems much more probable than a double photon or electron excitation. It is
a straightforward procedure to disentangle experimentally these two excitation mech-
anisms. In the course of the EELS experiment on C60 [6] multipole ionized states
up to C354 have been observed, with ionization states 1 to 4 in the ratio of 1 : 0.6 :
0.1 : 0.01. The proper treatment of such processes is beyond the framework of the
present paper and we plan to address these questions later. It will be highly desirable
to have experimental spectra for the emitted electrons, not only the ion yields as it

is presently the case, in order to be able to disentangle different possible channels.

 

Bibliography

[1] THE ALM’OST COIWPLETE BUCIXWIINSTERFULLERENE BIBLIOG-
RAPHY, available from F. A. Tinker, University of Arizona, Internet: tin—

ker@physicsarizonaedu .

['2] G. F. Bertsch, A. Bulgac, D. Tomanek and Y. Wang, Phys. Rev. Lett. 67,

26.90 (1.991).
[3] G. Barton and C. Eberlein, J. Chem. Phys. 95, 1512 (1991).

[4] I. V. Hertel, H. Steger, J. de Vries, B. Weiser, C. Menzel, B. Kamke and W.

Kamke, Phys. Rev. Lett. 68, 784 (1992).
[5] A. Bulgac and N. Ju, Phys. Rev. B 46, 4297 (1992).
[6] J. W. Keller and M. A. Coplan, Chem. Phys. Lett. 193, 89 (1992).
[7] N. Ju, A. Bulgac, and J. Keller, Accepted for publication in Phys. Rev. B.’

[8] J. H. Weaver, J. L. Martins, T. Komeda, Y. Chen, T. R. Ohmo, G. H. Kroll,
N. Troullier, R. E. Haufler and R. E. Smalley, Phys. Rev. Lett. 66, 1741

(1991).

[9] Y. Saito, H. Shinohara and A. Ohshita, Jpn. J. Appl. Phys. 30, L1068 (1991).

62

63

[10] G. Gensterblum, J. J. Pireaux, P. Thiry, R. Caudano, J. P. Vigneron, P.
Lambin and A. A. Lucas, Phys. Rev. Lett. 67, 2171 (1991).

[11] K. Zeppenfeld, Z. Phys. 243, 229 (1971).

[12] U. Diebold, A. Preisinger, P. Schattsdchneider and P. Varga, Surf. Sci. 197.
430 (1988).

[13] P. E. Batson, Phys. Rev. Lett. 70, 1822 (1993).

[14] D. Tomanek and M. Schliiter, Phys. Rev. Lett. 67, 2331 (1991).
[15] W. Kohn and L. J. Sham,iPhys. Rev. A 140. 1133 (1965).

[16] P. Hohenberg and W. Kohn, Phys. Rev. B 136, 864 (1964).

[17] D. Tomanek and M. A. Schliiter, Phys. Rev. Lett. 56, 1055 (1986).
[18] J.C. Slater and CF. Koster, Phys. Rev. 94, 1498 (1.954).

[1.9] S. Saito, and A. Oshiyama, Phys. Rev. Lett. 56, 1055 (1986).

[20] R. C. Haddon, L. E. Brus, and K. Raghavachari, Chem. Phys. Lett. 127.
459 (1986).

[21] M. Oziki, and A. Takahashi, Chem. Phys. Lett. 127, 242 (1986).
[22] S. Satapathy, Chem. Phys. Lett. 130, 545 (1986).

[23] P. W. Fowler, P. Lazzeretti, and R. Zanasi, Chem. Phys. Lett. 165, 79
(1990).

[24] L. D. Landau, Quantum Mechanics : Non—relativistic Theory, 3rd ed. (New

York, 1977).

 

64

[25] A. Bohr and B. Mottelson, Nuclear Structure, vol. II, Appendix 2, (Ben-

jamin, New York, 1974).

[26] W. A. de Heer, W. D. Knight, M. Y. Chou and M. L. Cohen, Sol. Stat.
Phys. 40, 93 (1987).

[27] S. Leach. M. Vervloet, A. Despres, E. Bréheret, J. P. Hare, T. J. Dennis, H.
W. Kroto, R. Taylor and D. R. M. Walton, J. Chem. Phys. in press.

[28] M. Braga, S. Larsson, A. Rosen and A. Volosov, Astrom. Astrophys. 245,

232 (1991).
[29] L. Goodwin, J. Phys. Cond. Matt. 3, 3869 (1991).

[30] C. Pan, M. P. Sampson, R. H. Hauge and J. L. Margrave, J. Chem. Phys.

95, 2944 (1992).

[31] S. Lundqvist and N. M. March, Theory of the inhomogeneous electron gas,
P. 180, ( Plenum Press, New York and London).

[32] H. Esbensen and G. F. Bertsch, Phys. Rev. Lett. 25, 2257 (1984).
[33] K. Yabana and G. F. Bertsch, private communication.
[34] K. Yabana and G. F. Bertsch, preprint (1993).

[35] M. Fieber—Erdmann, T. Drewello, W. Kratschmer and A. Ding, in Proc.
ISSPIC—6, September 15 — 22 (1992), Chicago.

 

Chapter 4

Finite Temperature Properties of
Small Sodium Clusters

In this chapter, we will present geometric and energetic properties of sodium clusters
with 8, 14, 20, 30 and 40 atoms at finite temperatures using an effective many—
body interaction among sodium atoms in the framework of an improved isothermal
molecular dynamics approach. We shall show that these clusters undergo two phase
transitions and that the two transition temperatures go up with the cluster size. The
two phase transitions are the equivalents of bulk melting and of boiling in a finite
system. At low temperatures, these sodium clusters are essentially incompressible,
but relatively easy to deform. At high temperatures they become extremely soft.
However, strong finite particle effects are observed. In particular, these clusters show

a more pronounced thermal expansion than the bulk.

4.1 Introductory Remarks

There are several reasons why it is interesting to study finite temperature properties
of atomic clusters, in particular of simple metal clusters. The most obvious one is the
fact that experimental results obtained so far are for relatively hot clusters [1, 2, 3].

Even though the experimentalists do not have at the moment a way of determining

65

66

the temperature of the clusters they produce, and only relatively rough estimates are
available, it is very likely that these clusters are melted. On the other hand, most
of the theoretical studies are performed at zero temperature. The properties and the
behaviour of finite systems at nonzero temperatures are of unquestionable theoretical
interest. While in the bulk there are well defined phase transitions (melting and boil-
ing) and statistical physics concepts (such as temperature, microcanonical, canonical
ensembles, etc.) are well defined, the validity of such terminology for finite systems
still raises some eyebrows, in spite of the fact that the usefulness and correctness of
such an approach have been justified for quite a while [4, 5, 6, 7]. Although the study
of finite systems opens an avenue towards understanding the way a piece of bulk ma-
terial emerges this field has its own set of problems. The relative role of the surface
versus volume effects is significantly enhanced compared to the bulk. An important
aspect is the presence of at least two types of degrees of freedom — electronic and
ionic — with vastly different characteristic time and correspondingly energy scales.
This fact seems to lend support to the validity of the traditional adiabatic approach.
Since the ion cores are so much heavier than the electrons and the physical phenom—
ena of interest occur in a region where the density of states for the ionic degrees of
freedom is extremely high [8], a classical description of these degrees of freedom seems
more than reasonable. One the other hand, a quantum description of the electronic
degrees of freedom is unquestionably the right approach. This View is supported by
the experimental evidence of pure quantum phenomena, reminiscent of the similar
nuclear phenomena: relative abundances of different atomic clusters, odd—even parti-
cle effects, geometrical shapes of the clusters (as indirectly inferred from the study of
their electromagnetic and optical properties), occurrence of the so called supershells;
as well as by simple theoretical estimates of the relevant electronic energy scales in

such finite objects [4, 5, 6, 7].

67

We attempt to elucidate to some extent the role of a finite temperature in sodium
clusters. The approach chosen by us is not irreproachable. In choosing one or another
description, various authors emphasize different dominant physical aspects. One can
distinguish roughly two dominant trends in the literature: 2') pure adiabatic picture
with an explicit treatment of the ionic degrees of freedom, which can be subdivided
in two subclasses ia) explicit treatment of the electronic degrees of freedom and 25)
implicit treatment of the electrons via an effective many—body classical potential for
the ions; and ii) featureless ionic background and explicit finite temperature treatment
of the electronic degrees of freedom. We shall not discuss here the quantum chemistry
calculations, which because of their rather ambitious goal to compute everything with

essentially no approximation, are limited to relatively small clusters [9].

The approach ia) seems to be accepted as the most “fundamental” way to describe
atomic clusters at zero as well as at finite temperatures and is often referred to as an
ab initio calculation [10, 11]. It amounts to a classical description of the ionic degrees
of freedom — which seems to be a very reasonable approximation, with the exception
of very light atoms like hydrogen, helium and to some extent lithium — and a Density
Functional Theory within the Local Density Approximation (DFT—LDA) for elec-
trons, augmented with a pseudopotential interaction between ionic cores and valence
electrons (which are the only ones usually treated explicitly) [10, 11, 12, 13]. From
the point of View of a pure practitioner, the main drawback of this type of approach is
the appalling amount of computer time required to extract relevant physical informa-
tion. It is for this reason that several groups decided to simplify the treatment of the
electronic degrees of freedom by using relatively old ideas, originated in condensed
matter physics, like the tight—binding approximation (see Refs. [15, 16, 17, 18] for a

sample of references), embedded atom [19, 20], etc. The “natural” extension of these

68

methods to finite temperatures still does raise some questions. Although one can
accept this type of approach for the description of insulators, where the gap between
the occupied and unoccupied electronic states is sufficiently large and the assumption
that the electrons are at all times essentially at zero temperature seems reasonable,
the applicability of the strict adiabatic approximation for metals or atomic metallic
clusters is questionable. One possible approach seems to be the generalization of the
DFT—LDA to finite temperatures [21]. However, this line of thought implies a rela-
tively quick relaxation of the electronic degrees of freedom. At low temperatures the
relaxation times in many—fermionic systems have a l/T2 behaviour [22, 23] and ionic
and electronic relaxation times can become comparable. Moreover, atomic clusters
are likely rather floppy objects and level crossing phenomena can occur. In such a
case nonadiabatic effects could prove to be quite important and in any event this

seems to be yet a term incognito.

A further simplification ib) relies on the use of an effective many—body interac-
tion among atoms, derived as before in the Born—Oppenheimer approximation for
the electrons [3, 24, 25, 26, 29]. In particular, the strongly delocalized character of
the valence electrons in simple metallic systems is at the origin of the many—body
character of this effective potential and it cannot be simulated by two—body forces
only. In this way one loses some pure quantum features, like shell effects, odd—even
effects, etc. However, one has the advantage, like in a Thomas—Fermi approximation,
of describing the gross properties of the clusters, since volume and surface effects

seem to be mimicked rather well by this type of model [25, 26, 27].

The last listed type of approach, ii ), emphasizes solely the role of electrons, while
the ionic degrees of freedom form a simple featureless jellium background, without

any dynamic or thermodynamic properties [1, 2, 3, 31, 33]. The main raison d’étre

69

for such models is the hypothesis that the electron dynamics dominates completely
the properties of simple metallic clusters, as manifested in abundances, odd—even ef-
fects, electromagnetic and optical properties, supershell effects, etc. While perfectly
suited to explain this range of phenomena, the jellium approach has its own inherent
limitations, especially when finite temperature properties of metallic clusters are con-
cerned. The gross properties of metallic clusters are rather poorly described in such
pure electronic models. For example the volume energy (total energy per particle)
derived in a jellium model is about 2.26 eV [31, 32] as contrasted with the cohesive
energy for sodium 1.113 eV. This fact alone sheds strong doubts on the treatment of
the ions as a physically featureless background. In the pure jellium model the shape
and stability of a cluster are totally governed by the electronic degrees of freedom.
If the volume energy is so strongly affected by the ions, a similar strong effect for
the surface tension, due to the ionic degrees of freedom, is plausible. As a result the
whole concept behind the jellium model, seemingly so successful in describing a whole

range of pure quantum effects, becomes to some extent questionable.

Moreover, a large series of molecular dynamics (MD) calculations [3, 5, 6, 7] point
to the existence of a large number of isomers, which are absent in a pure jellium
description, and their presence gives rise to the phase transitions occurring in atomic
clusters. The phase space for the ionic degrees of freedom increases exponentially
with increase of temperature (or equivalently excitation energy of the cluster, see
Section 3), while at the same time the phase space for the electrons increases at a
much smaller rate, since the electrons behave essentially like a degenerate Fermi gas.

One can estimate the electronic density of states using Bethe’s Fermi gas formula [34]

exp [(MZNE/sg (4 1)

pel(E): V \/Zl_8-E a

 

70

where N is the number of valence (delocalized) electrons in the cluster, 51: is their
Fermi energy and E is the electronic excitation energy. For the ionic degrees of
freedom an estimate from below (see Section 3) for their density of states can be
obtained from the classical formula for 3N —- 6 intrinsic degrees of freedom in the

harmonic approximation (Debye model) [8]
E3N— —7 3N

pion(E): (ZS—T—V—G)! 1:1} (433)
where N is the number of atoms in the cluster and w,- are the vibrational frequencies
of the normal modes. For the Nam cluster one obtains in this way p61 z 4 x 104
eV‘1 and pm, a: 1058 eV‘1 at an excitation energy of 3 eV. Consequently, the finite
temperature properties, the phase transitions in particular, should be dominated
by the ionic degrees of freedom. Since the density of states is directly linked with
the entropy, one can rephrase the above discussion in terms of entropic effects; the

electronic entropy is much smaller than the ionic entropy.

Experimentally, sodium clusters are likely to be produced in a liquid state, with
an estimated temperature around a few hundred degrees. This could serve as an
explanation of the spectacular success of the jellium model. Being in a liquid state,
the ions are rather mobile and therefore the cluster can be easily deformed. If the
contribution to the surface tension, originating from the ionic degrees of freedom,
is significantly smaller than the electronic part, one can expect that the electronic
degrees of freedom (through the quantum shell effects) would play the major role
in determining the shape and the stability of the cluster. The surface energy for
liquid bulk sodium (0 = [0.699 — 3.18 x 10‘4(T — Tmegtingfl eV) [35] seems to be in
qualitative agreement with the value extracted from jellium calculations, O'jeuz'um( T =
400 K) = 0.5918 eV [33] (The total surface energy is defined as 0N2/3, where N is

the number of atoms in the cluster). However, the surface tension for liquid bulk

71

sodium shows a relatively strong temperature dependence, while the jellium model
seems to underestimate it by at least one order of magnitude in the temperature
range T = 0...600 K, doJ-euium/dT % 2. ..8 X 10‘5 eV/K compared to 3.18 x 10‘4
eV/ K for liquid bulk sodium [35] (We have estimated the temperature dependence
of the surface energy in the jellium model from the results for the liquid drop energy
in Ref. [33]). This comparison seems to indicate that a combined treatment of the
ionic and electronic degrees of freedom is desirable in order to understand the finite
temperature properties of these clusters. The apparent agreement between the jellium
prediction and the bulk value for the surface tension might very easily prove to be an

accident.

MD calculations with effective many—body potentials definitely fail to describe
the electronic shell effects. The amplitude of the so called shell—correction terms
(computed in the spherical approximation, which overestimates them as a result,
probably by as much as a factor of two) never exceeds 1 - 3 eV for clusters with up
to a few hundred atoms and a little bit more for larger ones [3, 33]. When converted
into energy per atom, the magnitude of the shell—corrections is relatively small, about
0.025 eV or 300 K for N = 40 (note however that this is still about 20% of the
potential internal energy for Na40 around T = 500 K, see Section 3, and therefore a
quite sizable correction). For larger clusters the contribution from these pure quantum
effects becomes even smaller, smaller than, or at least of the same order of magnitude
as, the uncertainty of the estimated temperatures of the clusters. Nevertheless, the
rather fine details of the abundance distributions seem to correlate qualitatively well
with the effects predicted on the basis of the shell—corrections [33, 36]. One might
therefore conclude that, if one is not interested in rather fine details, a pure MD

approach to these metallic clusters might provide a quite reasonable description.

72

There is one more argument in favour of using such a relatively simplified ap-
proach. A reasonable model for the description of thermodynamic properties of small
metallic particles has been suggested by Kubo and Gorkov and Eliashberg (see Refs.
[37, 38] for reviews of this approach and other relevant references). The specific heat at
small temperatures for a many—fermion system, in an independent particle model. is
shown to be determined essentially by the fluctuation properties of the single—particle
spectra near the Fermi level. Assuming either a Poisson or GOE type of level spacing
distribution, one can show that the electronic contribution to the specific heat at low
temperatures, T < 25F/N, for an even number of electrons, is either Ce; % 5TN/2cp
or 061 % 30T2N2/452F respectively. Here T is the absolute temperature, 5;: z 3.2
eV is the Fermi energy and N is the number of electrons. Throughout the paper
we shall use energetic units for the temperature (k3 = 1) unless otherwise explicitly
mentioned. The mean level spacing is estimated as 5 F / N. Even though these values
represent averages over an ensemble of clusters and not values for individual clusters,
they should provide a good estimate of the magnitude of the electronic specific heat.
For temperatures of interest in sodium clusters, the ionic contribution to the specific
heat is dominant, since CW, 2 (3N — 6) >> Ce; (see Refs. [39]). At relatively high
temperatures the electronic contribution to the specific heat can become comparable
to the ionic part (For T 2 25p/N the specific heat for the electrons, in the leading
1/N approximation, is C31 2 27r2TN/65p, irrespective of the character of the fluctu-
ations [37, 38]). Consequently, the structural changes (phase transitions) are likely
to be caused by the ionic degrees of freedom for not too large and moderate sized
clusters. This last argument is another way of comparing the electronic and ionic

densities of states, see eqs. (4.1—4.2).

In spite of all these semi—quantitative or qualitative arguments which we have

73

presented in favor of dealing with mostly the ionic degrees of freedom in metallic
clusters at finite temperatures, we do not mean to imply that the electronic degrees
of freedom are unimportant. On the contrary, we believe that a combined treatment
of both ionic and electronic degrees of freedom is warranted. In Refs. [13, 14] such
a problem has been studied to a certain extent, however from those results one can-
not make a judgement concerning the relative role of different degrees of freedom.
The quite ambitious method chosen there (DFT—LDA in conjunction with the Carr—
Parinello method), aimed at a complete description of the system, did not allow the
authors to clearly disentangle the specific role played by ions and electrons. The
relatively delicate interplay between these two types of degrees of freedom is likely
to lead to interesting phenomena, some of them we have partially alluded to above
(part of which seem to be beyond the present formulation of DFT at finite temper—
atures). Even though in the present paper we shall focus by default on the role of
the ionic degrees of freedom, we plan to extend our studies in the near future to a

comprehensive treatment of all degrees of freedom in these clusters.

In the following section we shall formulate explicitly our approach. In particular,
we shall present an improved isothermal molecular dynamics scheme, based on a
previous development [46, 41, 42] of the Nose—Hoover method [43, 44]. Section 3 is

devoted to the presentation of our results, followed by a short summary.

4.2 Isothermal Dynamics

The delocalized character of the valence s—electrons makes the alkali clusters quite
different from noble gas clusters. In the case of argon clusters one need only introduce
an effective two—body interaction among the atoms, typically the Lennard—Jones

potential. Electron delocalization makes a two—body interaction among alkali atoms

74

physically unacceptable. Instead of treating the electrons explicitly we shall use a
phenomenological interaction: the many—body alloy potential [24, 25, 26, 28, 30].
This potential has a many—body attractive part and a two—body repulsive part. The
parameters of this potential were determined from ab initio calculations with no free
parameters [25, 26, 27]. A similar type of interaction has been used in the description

of other metal clusters as well [3, 28, 29, 30].

The Hamiltonian describing the properties of a sodium cluster used by us is

Ecoh X
1 ‘- CI/P)\/ Zb

 

N p2
H = Et..+v= —*—+
22m (

 

 

N N 7‘7“ q N ,,._.
exp [—2q (4 — )] — exp {—19 (A — 1)] , 4.3
g 1%; f0 PV Zb ,2; 7‘0 ( )
J¢i 17,4:
where the following values for the constants have been adopted ECO), = —1.113 eV

for the bulk cohesive energy, Zb : 10.4 for the effective coordination number in the
bulk, p = 9, q = 3 and r0 = 3.66 A for the nearest neighbor distance [25, 26, 27].
The second term in the rhs describes the attractive part of the interaction and has
a many—body character, due to the delocalized valence s—electrons. The last term

describes the short range repulsion between sodium atoms.

Since we are interested in the thermal behaviour of the clusters, we have to in-
troduce the coupling to a thermal bath. In a recent study of Na7_9 clusters [39], a
cubic coupling scheme suggested previously [46, 41, 42] was used. This scheme was an
efficient algorithm to achieve ergodicity and also to allow a relatively fast exploration
of the phase space [45] It works well for small clusters such as Na7_9, but the cubic
term introduced in the equations of motion makes the integration of the equations of
motion rather difficult at high temperatures and for large clusters, since it requires a

relatively small time step. For example, the choice of the time step ranged from 2.0

75

X 10‘15 to 0.15 x 10‘15 sec for different temperatures in the study of Na7_9 clusters
[39]. Another drawback of the cubic coupling scheme is that the optimal coupling
coefficients are different for different temperatures and clusters [39]. This situation
can be quite annoying, especially in the case of relatively large clusters. In this paper,
we propose a new coupling scheme, similar to the one proposed for the description of
a Brownian particle [46]. We will make the dependence of the coupling coefficients
on the temperature and cluster size explicit. Furthermore, we can use the same in-
tegration time step for all cluster sizes and temperatures. In all our simulations, we
used 1 X 10‘15 sec as the time step and the simulations were carried out for 106
steps. Consequently the system was evolved for l x 10'9 sec for all temperatures.
This makes the length of the simulation about 7 times longer than in the previous

study [39] at high temperatures.

Our new coupling scheme to a thermal bath at temperature T is described by the
following equations of motion for the coordinates :r,, momenta pr,- and pseudofriction

coefficients 53, Cr and 6,, (here for the :r—components only)

 

 

5,, = Br: (4.4)
'm
. 8V 080 3 P—xi 380 [92' 760 [33'
xi = — a _ — _ 17 if; _ _ 3—11' 4.5
p ()(L‘i NLO €x_ p0 IVLOC p0 a0 TVLQ6 p8 ( )

 

51::

iN=1pxi
_ 1 4.6
NmT 3) ( )

poLo

 

 

[305110 NmTpo NmT Npo

 

a:

N _ {V .
1. = T(~ z) (.7)

2i]: 1 pri_ 3 Z:iil Pit) (4.8)
POLO N mTPo N P6
Similar equations hold for the y— and z—components. In the equations above, e0 ~

ngwf) is a constant with the dimension of energy, of the order of the energy cor-

76

responding to the highest frequency-in the system (one can call it the Debye fre-
quency, tap, for a cluster). L0 is a constant with the dimension of length of order
~ 1 A, P0 = \/2m—T is the average thermal momentum at temperature T and (10 is
a dimensionless constant of order one. oz, 6, 7 are dimensionless constants, cluster

independent and their temperature dependence is given by

I mL2

where to is the smallest characteristic time scale of the system, i.e. to ~ 27r/wD. [n

 

relations. (4.9) the number of particles N appears explicitly because the amplitude
of fluctuations of the terms in parentheses in eqs. (4.6-4.8) is proportional to N.
All the constants, whose magnitude we have specified only by the order of magni-
tude, can be varied within reasonable limits, without critically affecting the quality of
the simulation. However, relatively large variations of these constants might require
significantly longer simulation times. This might not affect the ergodic properties of
the equations of motion, and therefore in theory will lead to correct results, provided

that the simulation is long enough.

A few clarifications, concerning the nature of this type of coupling to a thermo-
stat, are in order. A complex system is always characterized by a rather wide range
of characteristic frequencies (modes). In isothermal MD the thermostat is coupled
mostly to some modes. In a 3—dimensional system the density of modes is largest at
high frequency, like in the phonon density in a solid. Consequently, in order to have
the most effective coupling to a thermostat, the characteristic frequencies of the ther—
mostat should be comparable with the Debye frequency of the system under study.
An analysis of the above equations of motion shows that this is indeed the case. A

quick thermalization implies also that energy is exchanged at a reasonable rate among

77

all the modes of the system (i.e. relatively small apparent relaxation times). At the
same time the thermalization of the slowest modes of the system is achieved only if
the total time of the simulation is larger than the characteristic time of these slow
modes and the intrinsic relaxation times of the system. It is possible to devise an
isothermal MD in which coupling to all modes (both slow and fast) is achieved by a
slight generalization of the type of couplings studied so far [47]. One may wonder as
well, why we have introduced a relatively larger number of pseudofriction coefficients.
Our experience [45, 46] shows that by increasing the number of pseudofriction coef-
ficients one can ensure a more efficient rate of exploration of the system phase space
(smaller apparent relaxation times) and hence avoid problems with lack of ergodicity.
In particular one can achieve a better convergence in the case of phase transitions
when the so called critical slowing down phenomenon sets in [45]. Moreover, in the
case of isolated systems one wants to make sure that there are no conserved quanti-
ties (whose presence signals the absence of ergodicity), in particular neither the total
angular momentum nor its direction ia a constant of motion. The price one has to
pay, the increased number of differential equations to be solved, is insignificant, es-
pecially for large systems. At the same time, a quicker decorrelation time among the
generated phase space configurations (which is the essential element of any ensemble

average procedure) ensures a much better overall quality of the simulation.

Provided that the equations of motion generate ergodic trajectories, one can re—

place the phase space average by the time average, which is much simpler to compute

(A(p,9)) = 271;) / d3di3quxp [-H—(§.’—q)] 409.9)
= $151.10 ifamanmn) (4-10)

where p,q stand for the momenta and coordinates, Z (T) is the partition function and

78

A(p, q) is an arbitrary observable.

In this new coupling scheme, both the center of mass coordinate and momentum
of the cluster are not constants of motion. When we discuss the thermal properties
below, we will refer everything to the instanteneous center of mass of the cluster and
present only the intrinsic thermodynamic properties of the cluster. In order to avoid
evaporation at high T (around and above the “boiling” phase transition), we have
added a linear restoring force every time a particular interparticle distance r,,- > 3R,

where R = roNl/3 and N is the number of particles in the cluster.

4.3 Thermal Properties

We shall follow Refs. [39] to classify the cluster characteristics into two categories:
thermodynamic properties (internal energy, specific heat, density of states, phase
transitions) and geometric properties (shape, rms radius, momenta of inertia, relative
bond length). The energetics and geometry of the cluster are strongly correlated and
their analysis is extremely helpful in understanding the thermal behaviour of atomic

clusters.

4.3.1 Thermodynamic Properties

During each time step we have monitored the kinetic, potential, rotational, vibrational
and total energies of the cluster. The total kinetic energy carries no useful information
about a system in the canonical ensemble and we have used it merely as a check of the
quality of our simulation. However, the kinetic energy of a cluster can be separated
in an unique way into two nontrivial parts, rotational and vibrational energies [48].
Although one might expect that the rotational energy could be different from that

of a rigid body, since the clusters are to some extent floppy objects, in the whole

79

range of temperatures studied we have found that these clusters behave essentially
like rigid bodies. In particular the rotational specific heat of the clusters is CM, 2 3 / 2
as for a rigid body. This is due to the fact that in this temperature range the thermal

rotational motion is not fast enough to lead to any significant centrifugal stretching.

The only pertinent physical information comes from the analysis of the potential
energy of the clusters. For each temperature, we bin the value of the potential energy
at each time step, and construct the histogram f(E-p0t,T). This distribution of the

potential energy at a given temperature T is

 

1 Epot
f(E,,o,,T) = ——(—T)-p(Epot)exp (— ). (4.11)

Zpot T

where ,0(Epot) is the density of states originating form the potential energy only, and

 

ZPO,(T) = /,0(Epot) exp (- E51“) dEpot (4.12)

is the potential energy partition function. The distribution f(Ep0t, T) is peaked near
the average potential energy for the corresponding T, and drops off rapidly on both
sides. This allows for a reconstruction of the density of states over an energy range
comparable to the width of f(Epot,T) [39, 49]. By piecing together parts of p(Epot)
from simulations at different temperatures one can reconstruct the density of states
over a significant energy interval, up to an undetermined multiplicative factor, which
can in principle be determined as well. The logarithm Of the densities of states of the
potential energy for the clusters studied are shown in fig. 4.1 along with the logarithm

of the corresponding Debye density of states, eq. (4.2), divided by (3N — 8)/2.

One can clearly see an increased density of states at higher excitation energies,
which we associate with the softening of the clusters and which leads to the onset of

the phase transitions. The internal energy and the specific heat can be determined

 

80

 

 

 

    

VirVYfV'V"

N O o? GIN (a) aF U!
I .-.....

 

'V'VVVVVVI

   

ufi

"V'T'v

N

"V

 

Logic p(Epo.)/[(3N-6)/2-.1]

 

D

 

 

HNUaFCtho-F

 

 

Figure 4.1: The decimal logarithm of the unnormalized density of states (potential
energy only) as a function of the excitation energy per particle (solid line) and the
Debye model expectation (dashed line), divided by (3N - 8)/2.

81

from the following standard relations

 

 

1 E 0

U(T) :— <Epot> = m/Epotp(Epot)eXp (—‘ lit) dEpota (413)
< E2 > — < E o >‘2

C(T) = p“ N , (4.14)

T2
once the density of states is known. We have computed these quantities using the
above relations and checked them also against the corresponding averaged quantities
obtained during each separate simulation, and the agreement between the two meth-
ods was satisfactory. However, in spite of the fact that our simulation time is rather
long (1 nsec), the temperature dependence of some quantities is often not very smooth
and some fluctuations are still present (their magnitude will be evident in some of the
quantities presented in the following subsection, for which an equivalent of the density
of states cannot be defined). This is indicative of the fact that even longer simulation
times and/ or an optimization of the coupling scheme to the thermostat are needed.
The extracted density of states proved to be a rather smooth function of the energy
(the inherent statistical fluctuations, due to finite simulation times, are indistinguish-
able in the plots). The internal energies presented here are calculated with respect to
the ground state energy of each cluster. The temperature dependences of the internal
energy and of the specific heat are displayed in figs. 4.2 and 4.3. Even though the
change in the slope as a function of temperature is not so evident for U (T), one can
clearly identify the existence of two phase transitions in fig. 4.3. The lowest phase
transition can be associated with the melting and the highest one with the boiling of
the cluster. As expected, in finite systems the phase transitions are not sharp and not
as well defined as in an infinite system. However one can unmistakably conclude that
“something quite dramatic happened” to the cluster. As the cluster size increases,

the two phase transitions become better defined and in the thermodynamic limit they

 

 

 

 

 

82

will likely correspond to the melting and boiling phase transitions for bulk sodium,
with Tmeztmg = 371 K and Tbomng = 1156 K. The reason for the appearance of these
phase transitions is the dependence of the density of states on the excitation energy
of the cluster. At low excitation energies, below the melting point, the density of
states has a power law behaviour, characteristic for an ensemble of harmonic oscilla-
tors, see eq. (4.2). The cluster resembles a small crystal and the atoms perform only
small amplitude oscillations near their equilibrium positions. This is reminiscent of
the Einstein model for a solid. Note however that the finite number of ionic degrees
of freedom leads to a finite number of vibration energies. Once the excitation energy
is increased, the power law behaviour changes to an exponential one, see fig. 4.1,
the available phase space increases at a tremendous rate and as a result so does the
entropy of the cluster. At approximately 0.4le eV excitation energy, the density
of states for the cluster exceeds the Debye density of states, see eq. (4.2), by about
( 3N — 6) / 2 orders of magnitude. The logarithm of the density of states is essentially
the entropy of the system (5(E) = In p(E)). Our results shows that at this excitation
energy the entropy per particle exceeds the Debye model estimate by approximately
3-4 units. One can see that the two. transition temperatures go up with the cluster
size and both melting and boiling occur at lower temperatures in a cluster than in
the bulk [3, 5, 6, 7]. In finite systems these temperatures are obviously not sharply
defined and both melting and boiling occur in a range of temperatures. A cluster has
a significant part of its atoms (or perhaps better to say, most of them) at the surface.
These atoms are not strongly bound, since their mean coordination number is lower
than that for the inner atoms, and it is relatively easy to “shake them up”. For that
reason atomic clusters “below the melting point” are in a rather unusual state, some-
times called the glassy, the molten or the fluctuating state [3, 28, 29, 30, 39]. “Above

the melting point” they start behaving like liquids, as one can see more easily in the

 

83

behaviour of their pair correlation function (see the next subsection and Refs. [39].

As a side remark, we would like to stress the advantage of performing an MD
calculation in the canonical rather than microcanonical ensemble. One finds often
in the literature the statement that the two approaches are essentially equivalent
with respect to the amount of physical information extracted and for that reason
different authors rather often prefer microcanonical simulations to canonical ones. In
a canonical simulation one can relatively easily extract the density of states, which
is unaccessible in a microcanonical one, and the thermodynamic behaviour of the
system can be easily inferred and understood. For a finite system. the density of
states plays a similar central role as in statistical physics of large or infinite systems.
In particular, it is much less expensive to perform a canonical MD and find structural
isomers than in a straightforward search. If an isomer exists, its presence will show
up at sufficiently high temperatures, as, in particular, a more careful analysis of the
behaviour of the density of states shows (To some extent a similar ideology is behind

the popular simulated annealing method).

In spite of the fact that we are dealing here with such small systems, we decided
to estimate the latent heat of fusion for sodium from these simulations. Since the
first phase transition occurs around 200—300 K, and the transition peak is very broad,
we will simply take 500 K as the temperature at which we regard the clusters are
completely in a liquid state. If there is no phase transition, the potential internal
energy at temperature T would be (3N — 6)T/2, as for pure harmonic oscillators
(note that the ionic degrees of freedom are treated classically). We shall identify the
excess over the internal energy of pure oscillators as the energy necessary to melt
the cluster, i.e. the energy needed to cause a structural change of the system. We

estimate the heat of fusion to a value between 2.36 and 3.45 kJ / mole, which compares

S~l

 

 

 

-A LALAAJAA

l'

LAJA- A

 

 

-1 l m..

1000 1250

 

 

o- --IAALmlLAAAI-
0 250 500 750

T [K]

Figure 4.2: The temperature dependences of internal potential energies per particle.

S 7)

 

 

 

<111<<1<1141d111 11111144 111111441114111111 414411111111441111141 11d 11.4141444141 1411114441411111
1 l 1 1 I.
A A A A A
v A A A A
v A A A A A
A A A A A
.I I IA I. L I.
. A A A A A
. A A A A A
v A A A A A
. A A A A A
I 1 1 IA IA 1
v A A A A
v A A
v A A A A A
r A AA A A A
T l l 1 IA
4 A A A A A
v A A A A
. w A w a. A m A 9“ A a A
. a A a A a A a A a A
.I N 1 N IA N L N IA N A
. A A . A A A
v A A A A A
A A A A .
A . . A A
butt»p-ypu-nbpnp-npnp-Dpp nthbbnbhhbbbhbbrbh ppbbhpbbbbbbbbbbbbbbb bbbbbEbbbhbbe- rp-hpbtnbbbb-bPD

 

 

 

 

 

505
221.

750 1000 1250
T [K]

500

250

Figure 4.3: The temperature dependences of speclic heat.

86

unexpectedly well with the value 2.6 kJ / mole for bulk sodium [35]. One can try to
estimate in the same manner the heat of vaporization as well. At high temperatures
the quality of our simulations is worse, due to significant evaporation. Nevertheless,
we obtain a value which is approximately equal to the experimentally measured heat
of vaporization for sodium 89.6 kJ/mol (within 15 — 20 ‘76). This is not very surprising,
since the effective interaction we are using has the correct asymptotic behaviour of
the cohesive energy for large systems built in, and the heat of vaporization is almost

equal to it.

4.3.2 Geometric Properties

A very sensitive indicator of structural changes in a cluster is the rms relative bond

length [5, 6, 7, 39]

 

 

2 N <r,-2->—<r,-->2
éz—Z“ ’ ’ . (4.15)
1V(1’V _1)i<j < Tij >

The temperature dependence of the rms relative bond length is shown in fig. 4.4. As
we have mentioned above, at low temperatures, “below the first phase transition”, the
cluster behaves like a small crystal, with the atoms oscillating with small amplitudes
around their equilibrium positions. At temperatures around 200 K the atoms are at
their equilibrium positions for relatively long periods of time. However from time
to time an atom jumps from one equilibrium position to another clue to thermal
fluctuations [3, 39]. The rms relative bond length changes drastically around this
temperature, indicative of a structural change. At still higher temperatures, the

atoms become extremely mobile and move across the entire cluster [39].

In the region of the second phase transition 850—1000 K, we observe a second

increase in 6. This is partially linked with the fact that some atoms can evaporate,

 

 

 

 

 

 

 

 

 

    

 

 

0.635Ti'h PF fx ' .1
0 5- Na x x x '
:35. x x x x x i
0.2E-Ax l l AlAAA 1AAAAL‘:
0.8} ‘ X X x ‘1
I X .
0.5? N330 x «I
0.43- x l
03E x x x x .5
E x ' 3
0.2:A AAA AALAAAALAAAA;AAAAA
0.6? Nam x X X X '3
0.5? x ‘i
'o E i
0.4:- x X 1
cat X X X <1
. E x l A l A IA AlAA All
i X x 1

. :- x x ‘2
06E N314 x :
0.4; x x 1
0.3;- x x x «;
t X 1
0.2 A 1AAAAAAA AAA AAAAA 1‘
x x x .

0.9 Na. " x 1
x 1

0.4 1
x X X X X ii

0.2 A¥AALAAAA1AAA LAAAAIAAAALl
0 260 500 760 1000 1250

’1‘ [K]

Figure 4.4: The temperature dependences of the relative bond length 6.

 

88

even though at a not very significant rate yet. At temperatures above 1000 K, the
evaporation of atoms is significant and the precision of our results is therefore re-
duced. A better approach will be a grand canonical ensemble at these temperatures,
suited for the description of the liquid—vapor coexistence. The structural transition
in these sodium clusters can be seen also in the behaviour of the (unnormalized) pair
correlation function f(r,-j), see figs. 4.5 and 4.6. Below the first phase transition
the interparticle separations are rather well defined, as one would expect for a small
crystallite. In the fluid phase, even though one can see a well defined short range
correlation among particles, the long correlation is almost completely lost and at very

high temperatures the presence of the vapor phase is obvious.

For each spatial configuration of the cluster, we have calculated three geometric
quantities: the rms radius of the cluster, which characterizes the cluster size ; the
shape parameter 6, for the degree of asphericity of the cluster and the shape parameter
’7 (0 S 7 g 71' / 3), which describes the triaxiality of the cluster [39]. These parameters
provide an average information about the size and shape of the clusters. The rms
radius and shape parameters 6 and 7 are related to the principal momenta of inertia

(11 Z [2 Z 13 Z 0) through the following relations

k_
1,,=§7~’2 [1+fisin (7+L46—3M)] k:- 1,2,3, (4.16)

where

1 2 _
r: NZ!” (Zn-.0) (4.17)

is the rms radius of the matter distribution. The condition that the semiminor axis
of the associated ellipsoid of the inertia is positive

2
2 T‘

a:—
3

[1 — 2fisin(7 + 7r/6)] 2 O (4.18)

f (rij) (arb. unit)

1.0

0.5

0.0
1.0

0.5

0.0
1.0

0.5

0.0

89

 

 

 

 

 

  

 

 

 

 

 

 

b I ' T f I ' ' I V f ' f l ' f ' i l ‘

l

D

D

’ T=1000 K N 20 ‘

- a «

D

' 1

D

’ 1

" -1

’ 1

i 1

' 1

’ 1

" :1

D

D

l 1

D

. -1

D

D 1

D 1

l A A A1 A A A A AAA- 1 A A A A 1 A L 1

D

D 1

D ' 1

D — 1

,. _ a ,

D 1

D 1

D

’ 1

l' ‘4

D 4

D

D 1

' 1

b '1

l l

D

D

D 1

D d

p

D

i l

l A A A 1 AA A A A A A A A 1

I. I U T V ' .1

D “

D

D 1

’ T-lOO K N ‘

D - a q

k 1

D

D

L 1
<

’ 1

D 1

' 1

D 1

D d

D

’ 1

D 1

D 1

p ‘1

D 1

D

' 1

’A A A A 1 A A AL 1 ‘

 

rij [Aug]

Figure 4.5: The pair correlation function for Nam at several temperatures.

f (rij) (arb. unit)

1.0

0.5

0.0

1.0

0.5

0.0
1.0

0.5

0.0

90

 

 

 

 

   

 

 

 

 

 

 

 

 

Figure 4.6: The pair correlation function for Nam at several temperatures.

P f I ' I - - l ' v V ' I f f ' T '

D
D
D
’ 1:1000 K N 40
" a
D
D
D
D
p
D
D
D
D
D
D
D
D
D
D
D
D
f A A A A A A L A A ‘A A A A L L
p ' V r I’ r T I v v ' v ' f I f ' T .
D
* T-5OO K/Na40
. -
D
D
D
D
b
D
D
D
D
D
b -
D
D
D
D A A k L A l
P ' T v v f V ‘0 v v f T— v ‘
f 1
f 1
i 1
f f
i 1
f 1
D
D
D -1
f l
f 4
i 1
f 1
p '4
f 1
f 1
D
i 1
- cl
[ 1
D
D

A 1A A A l A L

91

defines the region of allowed values for [3 and 7. If 7 = 0, then 11 = 12 > [3, and the
shape of the cluster is an axially symmetric prolate ellipsoid (cigar) with O S ,8 S 1.
For 7 = 7r/3, 11 > 12 = [3 and the shape corresponds to an axially symmetric oblate
ellipsoid (pancake) and 0 S 5 S 1/2. For all the remaining angles in between, the
shape is a triaxial ellipsoid. For fl = 1/2 and 7 = 7r/3, the shape is a disk of zero
thickness, while fl = 1 and 7 = 0 corresponds to a linear chain. The region of allowed
values for [3 and 7 has the shape of a triangle in the plane where [3 is the radial

distance and 7 is the angle at origin, see fig. 4.7.

The temperature dependence of the average principal momenta of inertia and
its relative covariances are displayed in figs. 4.8—4.9, the rms radius of the cluster
1", and the shape parameters 5, 7 along with their corresponding covariances, are
displayed in figs. 4.10—4.12. From these plots, one can see that all these clusters
tend to acquire a cigar—like shape at high temperatures. This might seem a bit
peculiar, but is quite easy to understand. One can show that the shape space measure
is proportional to i?“ sin 37(1/3d7 [52]. Consequently shapes with small values of 1'3
are strongly suppressed and for the same reason pure oblate or prolate shapes are
suppressed as well. At the same time the shapes with as large values as possible for
B are strongly favored by this measure and for this reason prolate—like shapes (for
which large values of fl, 1/2 S B S 1, 0 S 7 S 7r/6 are possible) start dominating
at high temperatures. For 0 S 3 S 1/2 oblate—like (7r/6 S 7 S 7r/3) and prolate—like
(0 S 7 S 7r/ 6) shapes are equally probable, according to the shape measure only.
One has to keep in mind that these distributions reflect the shape of the free energy
of the cluster in the shape space, which takes into account both energetic and entropy
properties of the system (the shape space measure is already included). We could have

defolded the shape measure from these plots (this is the standard convention in nuclear

 

 

 

 

Figure 4.7: The allowed region for the shape variables 13 and 7.

 

 

 

 

 

 

 

 

93
150: 'l I r 1

125:- °;
100} Nam o x;

: O X
75? X + 1
50: 2 Q9 + + “
25issxx޴+++ .
150EAAATLT¥ l -hl 1 L:
1255- i
: 03
100;- Naao o x«j
75? fix 1
50;- 29 +4. -
25'n!¥'¥¥3¥++ i
F-'100 A1 mil 1 in L61:
N 80:- o 1
° : Nazo O X?
360;- x x ~;
:40? g2 + +';
22- + ~‘
_°5¥¥$88¥§++ 3
lootk 1 h; 1krgrrr- -1 - I;
80? Nal‘ 0'?
60;" ng'l
40? 2X
205- +4
Enum$¥8$++++i
60’: 1 --A l A 01:
50;- o -f
40_ N33 2Xx;
3°? 2 '
ZOE" 2 1
:- +1
135A:glgmgh¥A$A¥1f-f--r-f*‘li
0 250 500 750 1000 1250

T [K]

Figure 4.8: The temperature dependences of the average principal momenta of inertia
per particle.

 

 

 

 

 

 

 

 

 

 

m
0.5,”-.fiver -, :
0.4; N ‘ 8 + '3
: a +
0.3:- 4° X 5 + 5 z, .
E o +
0.2;- x + ~
0.1% . u 2 .
00:13 !W“-
ole:- .
I N330 5 3 3 + i
0.4} + 5 +.
~ E + " i
0.2.,b g 7
ooEl.-*l!-¥ L? l k-l - Li
'I + 1
0.6? _ x + 1
I N320 g 3 + 1
H04; - ‘ § '1
E '5 x . 3
0.2:- § <1
I I F 1
. n 1‘ .
0.04 L-h--1-- -1----1- 1‘
0'8? X + + i
0.63- Na“ 2 3‘ +4
0.4:- 3 l' 24
0.23- «I
ooLL‘Lg-¥-¥-¥-kL-Lk 1--- 1
120i- x + + .
08:. N36 0 6 ‘l' ..
0.3;- + , -
0.4? E + *1
8.31:4 -¥¥-¥-¥--1A---1-- 1‘;
' o 250 500 750 1000 1250
TlK]

Figure 4.9: The temperature. dependences of the relative covariances of the average
principal momenta of inertia. '

95

physics calculations of deformed nuclei [3, 52], but since any average will involve this
measure anyway, such a way of displaying the shape properties of the cluster seems
to us to a certain extent misleading. At low temperatures, these clusters are almost
incompressible (their rms radii have relatively sharp distributions), but at the same
time they can change their shape rather easily (the distributions for both fl and 7
have significant widths). At high temperatures, these clusters become apparently
soft as well, having the characteristics of a compressible fluid. This apparent softness
is to some extent an indirect indication of the presence of the vapor. In spite of
the shape differences of their ground geometries, the high temperature behavior of
their shapes is quite similar. Thelcharacter of the shape parameters distributions
changes dramatically in the Vicinity of the phase transitions. The smallest principal
momentum of inertia has the biggest fluctuations, since atoms evaporate more readily

from the sides of the cigar for obvious reasons (more chances).

We have tried to extract the coefficient of thermal expansion from our results.
Unfortunately we did not have enough statistics for a precise determination of the
linear and quadratic expansion coefficients. Nevertheless, as one see from fig. 4.10
the clusters display a rather well defined nonlinear thermal expansion. This fact
might play a quite significant role in explaining the red shift of the Mie plasmon in
sodium clusters [1, 2]. In the jellium calculations, routinely used in describing the
optical response of sodium clusters, the bulk density for the jellium is assumed. Be-
sides melting and boiling in a quite different way than the bulk, sodium clusters seem
to have also a rather distinct thermal expansion behaviour. Since the experiments
are lilcely performed with liquid clusters, their larger volume can easily serve as an
argument in favor of a lower plasmon energy. The classical Mie plasmon frequency

will change approximately as 6wM,e/wM,e = —36r/2r when the radius of the cluster

 

‘ 1: I: :1 q: :4 <44<r4r4u<44<q4<<<qr144 rd¢rrrd<441qr11q44¢4q4 <44:r1«<<<1<<44r<‘d<4 1r¢dr4<<d41¢ddr11q

96

New

 

 

r
AW
3

750 ‘
T [K]

L
k

500

{m m
g
.N N

xxx3‘

250

A

x.

x
X
X
m. x r
N
X
x

X
x A

 

 

 

 

4
r

. . A

 

14;-
12,
1o:-

 

20864208640864 aa‘z
11‘ 1.1 1

flag; .—

 

Figure 4.10: The temperature dependences of the rms radius and its covariance (error

bars). '

0.6 ’

0.6
0.4
0.2

0.0 ’

0.6
0.6
0.4
0.2

0.0 I

0.6
0.6
0.4
0.2

0.0 I

0.6
0.6
0.4
0.2

0.0
0.6

0.6
0.4
0.2
0.0

97

 

 

 

p f l I l I ‘
D
" 1
' Na“, ‘
’ 1
' 1
’ d
- 1
’ 1
h 1
' 1
’ d
- 1
’ 1
I 1
' 1
- 1
: 3 8 :
D

A A A A l A A A A l A A A A l A A A 1 A A A A l ‘
’ 1
’ 1
’ 1
. 1
’ Naao ‘
' 1
’ 1
’ d
b 1
’ 1
’ 1
. 1
’ d
- 1
’ 1
’ 1
' 1
- 1
I . ‘
b . ‘

x A A l A A A A l A A A A l A A A 1 A A A A l ‘
1

. 1
’ 1
D ‘ ‘
' I; 1
' Elm 1
’ 1
’ d
. 1
. 1
’ 1
’ 1
' fl
- I
I

xx$$$%

AlAAAAl AAA!

 

Na“ .
MHHHH

 

v“ vvv

vvv'

'vvv

 

x¥:i¥f%*¥%% ?

A. ‘

d
1

 

. . A

 

260600 7 1000 1260

'1‘ [K]

Figure 4.11: The temperature (lependem‘es of the shape parameter :3 and Its COVELI'I-

ance (error bars).

Figure 4.12: The temperature dependences of the shape parameter '7' and its covari-

ance (error bars).

1 [Deg]

 

 

 

 

 

 

 

 

 

 

 

 

 

50;- Na“0
4o;- ‘3
so;-
10g-
DELL LIAAAAI A 1A AAlA AAl:
403- N330 ‘3
30:-
20? Ԥ
10;- ‘Q
0; AA 1AA'AAL ‘AAArAAAArAAAAl:
301- - Nam 11
20} 1*
10;- i} A ‘
o; ‘ ‘ 41
40:-
so;-
zo;- :
10? :
o: ALAAA :
30%- :
20:: ‘
10;} ‘ ;
01-“ ,A WA .‘1
o 250 600 750 10 1250

1' [K]

 

 

 

99

changes. At 400—600 K (the estimated temperature of the clusters in nozzle exper-
iments) the linear dimensions of a cluster seem to be about 6 — 15% larger than at
T=0 K (and significantly larger for higher temperatures), already of the right order
of magnitude to explain the observed red shift of the Mie plasmon. At the same time,
at these temperatures the volume distribution of the clusters have a quite sizeable
width Ar/r % 0.04 - 0.07 (and significantly larger at higher temperatures), which can
explain at least part of the width of the plasmon. In figs. 4.13—4.13 we present the
spherical part of the ionic density, extracted from our simulations at several selected
temperatures (The relatively large fluctuations of p(r) near the origin are due to ob-
vious volume effects, near 7‘ z 0 the probability to find an ion at a distance r from
the center of mass of the cluster is proportional to -r2p(r) and correspondingly the
statistics is lower). With increasing temperature this density becomes flatter. The
radius of the cluster does not seem to vary significantly, however the surface diffuse-
ness is greatly increased after the first phase transition, which explains the thermal
behaviour of the rms radius discussed above. This relatively large apparent surface
diffuseness is partially due to the fact that the cluster is rather strongly deformed as
well. One can conclude that the anomalous thermal expansion of sodium clusters,

when compared to the bulk, is mainly due to these surface effects.

An additional contribution to the width of the Mie plasmon in sodium clusters
at finite temperatures will arise from shape fluctuations. The rough estimates we
presented here should be looked upon only as order of magnitude effects. The shape
of the cluster in its ground state and that with one Mie plasmon excited can be quite
different and the ultimate shape and width of the plasmon should be estimated using

a more detailed procedure [53].

 

p(r) (unnormalized)

100

 

v v v '

T=1000 K/NaZO

   

F

v—

T=500 K/Nazo

db

cu-
(
1
1
q

 

ffi r v

 

 

L g 1
' I

T=300 K/NaZO

A L
‘

AAI A L
Vfirw ff

T=100 K/Na20

‘MAA LLA

0.0A A 2.5 5.0 17.5
1' [Ang]

I.

  

 

10.0 12.5

Figure 4.13: The extracted spherical part of the radial (unnormalized) ion density

distribution for Nago.

 

15.0

 

p(r) (unnormalized)

lOl

 

 

 

 

 

 

 

 

 

 

 

 

 

T=1000 K/Na40

E
T=500 K/Na40

T=300 K/Na40 ‘:

l ---;+-ve-.-n-i

,L' ' - g

&_ T=100 K/Na40 4

0.0 2.5 5.0 7.5 10.0 12.5 15.0

r [Ang]

Figure 4.14: The extracted spherical part of the radial (unnormalized) ion density
distribution for Na“). '

102

A comment concerning the potential use of such distributions in jellium—type
calculations is in order. These distributions represent ensemble averages. The time
scales for the ionic degrees of freedom are significantly lower than for the electronic
degrees of freedom. If one would like to include such thermal effects in the calculations
of the electronic properties and optical response of a cluster, one should perform the
corresponding jellium calculations for one member of the ensemble at a time and take
the ensemble average only afterWards. In particular, spill—out of the electron density

can be more pronounced for some members of the ensemble.

4.4 Closing Remarks

Since the energy scale characteristic for the ionic degrees of freedom is relatively small,
a huge phase space becomes available upon increasing the temperature of the cluster.
These entropic effects show up in rather drastic structural changes and, without being
extremely pedantic concerning the adequacy of the terminology, one can characterize
these changes as phase transitions in finite systems. With increasing temperature,
the clusters go first through a glassy/molten or fluctuating state (200—300 K) and
eventually become totally liquid (above 300 K). In this temperature range they are
almost incompressible, but highly deformable. At still higher temperatures (850—1000
K) they start boiling, vaporization sets in and they become rather soft / compressible.
In contradistinction with noble gas clusters [5, 6, 7], the melting and boiling do not
seem to be a geometric or particle number effect and the same behaviour is observed
for all the clusters we have studied (N = 7, 8, 9, 14, 20, 30 and 40) [39]. The thermal
expansion properties of these clusters seem to be more pronounced when compared to
the bulk and it is mostly due to an increased surface diffuseness and deformation of the

cluster. Thermally induced rotation is never fast enough to lead to any rotationally

103

induced effects and in the temperature range studied, one can safely characterize
these sodium clusters-as rigid bodies with respect to rotational degrees of freedom.
The onset of the new phases can be. observed in both thermodynamic and geometric

properties of the clusters, which have a strong temperature dependence.

The present approach does not account for the electronic degrees of freedom ex-
plicitly and therefore electronic shell effects are not accounted for. We expect however
that, due to the large ionic entropy effects we observe, similar behaviour should be
observed in a more complete description of these clusters. The electronic shell effects
should be rather important at relatively low temperatures, below melting. It will be
extremely interesting to study the effects of the large geometric fluctuations on the
electronic properties and on the optical response in particular (unless in a combined
treatment the clusters do not become more stiff. which does not quite seem to happen

[1:3, 14]).

Bibliography

[1] W. A. de Heer, W.D. Knight, M.Y. Chou and M. L. Cohen, Sol. Stat. Phys.
40,93 (1987).

['2] W. A. de Heer, to be published in Rev. Mod. Phys. (1993).
[3] S. Sugano, rWicrocluster Physics, Springer—Verlag (1991).

[4] T. L. Hill, Thermodynamics of small systems, parts 1- — 2, W.A.Benjamin.
New York and Amsterdam, (1963 — 1964).

[5] R. S. Berry, T. L. Beck, H. L. Davis and J. Jellinek, Adv. Chem. Phys. 70,
75 (1988).

[6] H. L. Davis, J. Jellinek and R. S. Berry, J. Chem. Phys, 86, 6456 (1987).
[7] R. S. Berry, J. Jellinek and G. Natanson, Phys. Rev. A 30, 919 (1984).

[8] G. F. Bertsch and R. A. Broglia, Oscillations in Finite Quantum Systems,

ch. 10, Cambridge University Press, in press (1993).

[9] V. Bonacié—Koutecky, P. Fantucci and J. Koutecky, Chem. Rev. 91, 1035

(1991).

[10] R. M. Dreizler and E. K. U. Gross, Density functional theory, Springer—
Verlag (1990).

104

105

[11] M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias and J. D. Joannopoulos,
Rev. Mod. Phys. 64, 1045 (1992).

[12] I. Moullet, J. L. Martins, F. Reuse and J. Buttet, Phys. Rev. Lett. 65, 476

(1990)

[13] P. Ballone, W. Andreoni, R. Car and M. Parrinello, Europhys. Lett. 8, 73

(1989)
[14] U. Réthlisberger and W. Andreoni, J. Chem. Phys. 94, 8129 (1991).

[15] D. J. Chadi and R. M. Martin, Solid State Commun. 19, 643 (1976).

 

[16] D. Tomanek and M. H. Schliiter, Phys. Rev. Lett. 61. 2331 (1991).
[17] L. Goodwin, J. Phys. Cond. Matt. 3, 3869 (1991).

[18] R. Poteau and F. Spiegelmann, Phys. Rev. B 45, 1878 (1992).

[19] .J. K. Norskov and N. D. Lang, Phys. Rev. B 21,2131 (1980).

[20] .1. K. Norskov, Phys. Rev. B 26,2875 (1982).

[21] U. Gupta and A. K. Rajagopal, Phys. Rep. 87, 259 (1982).

[22] A. A. Abrikosov, L. P. Gorkov, I. E. Dzyaloshinski, Methods of quantum field

theory in statistical physics, Dover Pubs. (1975);

[23] D. Pines and P. Nozieres, The theory of quantum liquids, W. A. Benjamin
(1966). Note, however, that the l/T2 behaviour of the fermionic relaxation
time is strictly valid for infinite systems and its extrapolation to a finite

system should be taken with certain care.

[24] R. P. Gupta, Phys. Rev. B 23, 6265 (1981).

106

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

[26] D. Tomanek, S. Mukherjee and K. H. Bennenmann, Phys. Rev. B 28, 665
(1983); ibid 29, 1076, (1984).

[27] D. Tomanek, A. A. Aligia and C. A. Balseiro, Phys. Rev. B 32, 5051 (1985).
[28] S. Sawada and S. Sugano, Z. Phys. D 14, 247 (1989).
[29] J. Jellinek and I. L. Garzon, Z. Phys. D, 20, 239 (1991).

[30] I. L. Garzon and J. Jellinek, Z. Phys. D. 20, 235 (1991).

 

[31] w. Ekardt, Phys. Rev. B 29, 1558 (1984); and 31,6360 (1985).

[32] D. E. Beck, Phys. Rev. B 30, 6935 (1984); and 35, 7325 (1985).

[33] O. Genzken and M. Brack, Phys. Rev. Lett. 67, 3286 (1991).

[34] C. A. Engelbrecht and J. R. Engelbrecht, Ann. of Phys. 207, 1 (1991).

[35] “CRC Handbook of Chemistry and Physics”, 70th Edition, CRC Press, Boca

Raton, Florida, 1990.

[36] J. S. Pedersen, S. Bjornholm, J. Borggreen, K. Hansen, T. P. Martin and H.

D. Rasmussen, Nature 191, 733 (1991).
[37] W. P. Halperin, Rev. Mod. Phys. 58, 533 (1986).

[38] T. A. Brody, J. Flores, J. B. French, P. A. Mello, A. Pandey and S. S. M.
Wong, Rev. Mod. Phys. 53, 385 (1981).

[39] A. Bulgac and D. Kusnezov, Phys. Rev. Lett. 68, 1335 (1992); Phys. Rev.

B 45, 1988 (1992).

107

[40] A. Bulgac and D. Kusnezov, Phys. Rev. A 42, 5045 (1990).
[41] D. Kusnezov, A. Bulgac and W. Bauer, Ann. of Phys. 204, 155 (1990).
[42] D. Kusnezov and A. Bulgac, Ann. of Phys. 214, 180 (1992).

[43] S. Nosé, J. Chem. Phys. 81, 511 (1984); Mol. Phys. 52, 255 (1984); Prog.

Theor. Phys. Supp. 103, 1 (1991).

[44] W. G. Hoover, Phys. Rev. A 31, 1695 (1985); Molecular Dynamics, Lecture
Notes in Physics, Springer—Verlag, New York, (1986), Vol. 258, Computa-

tional Statistical Mechanics, Elsevier, Amsterdam (1991).

[45] J. H. Sloan, D. Kusnezov and A. Bulgac, Computational Quantum Physics,
AIP, New York, vol. 260, 23 (1992); Nucl. Phys. B 26, 623 (1992).

[46] A. Bulgac and D. Kusnezov, Phys. Lett. A 151, 122 (1990).

[47] A. Bulgac and D. Kusnezov, to be published.

[48] J. Jellinek and D. H. Li, Phys. Rev. Lett. 62, 241 (1989).

[49] P. Labastie and R.L. Whetten, Phys. Rev. Lett. 65, 1567 (1990).
[50] S. Iijima and T. Ichihashi, Phys. Rev. Lett. 56, 616 (1986).

[51] P. M. Ajayan and L. D. Marks, Phys. Rev. Lett. 60, 585 (1988); ibid 63,
279 (1989).

[52] A. Bohr and B. Mottelson, Nuclear Structure, vol. II, Appendix 2, Benjamin,

New York, (1974).

[53] G. F. Bertsch and D. Tomanek, Phys. Rev. B 40, 2749 (1989);

Chapter 5

Summary and Perspectives

5.1 Summary

Cluster physics is an exciting and young field, it has not only brought methods and
concepts from different disciplines into new and unanticipated combinations, but also
led to the development of new and advanced experimental techniques. Besides a few
words about our work on the nuclear curvature energy problem, we have devoted most
of chapter 1 to the review of the new developments in cluster physics, with special

attention paid to the theoretical side.

Chapter 2 has been devoted to the calculation of the correlation energy of spherical
nuclei using RPA method. We found that the correlation energy has a very different
A-dependent behavior from that of the mean-field ground state energy. By finishing
this project I learned the theoretical tools used for the study presented in chapter 3
even though by that time I did not know I would jump on the bandwagon of doing

cluster physics research.

Chapters 3 and 4 are the main part of this thesis. We calculated the multipolar
plasmon excitations in fullerenes in chapter 3. These calculations consist of two parts.

One is concerned with the plasmons in C60, its treatment, due to its spherical sym-

108

109

metry, is simpler than that of the nonspherical clusters. The other one is about the
plasmons in other nonspherical fullerenes, their treatment is somewhat more compli-
cated. With some manipulation of the RPA equation, we are able to completely solve
the RPA problem within the framework of tight—binding description for the single
particle spectrum and wave functions. The equivalents of both 71' and a plasmons in
graphite are predicted for these clusters. Good agreement is obtained between our

calculations and available experimental data.

We addressed the finite temperature effects of clusters in chapter 4. Since clusters
are produced at high temperatures, the study of clusters at finite temperatures is
very important and absolutely warranted. We have computed rotational, vibrational
and potential energies and corresponding specific heats, momenta of inertia, shapes
parameters and bond lengths of Na8,14,20,30,40 over a wide range of temperatures.
The equivalents of both “melting” and “boiling” of the bulk are predicted for these

clusters.

5.2 Perspectives

There are a number of directions along which the works presented in this thesis can
be extended. First of all, the extremely large width of the a plasmon in fullerenes is
an important problem to be understood. Several possibilities come to mind. As we
have remarked in Chapter 3, because 2t.)7r and 3t.)7r are very close in energy to 3w.”
there may be a strong, essentially resonant coupling between them. This coupling of
electronic excitations may lead to the large width of the a plasmons. Another route
one may take to improve the quality of the calculations in chapter 3 is to use a better
Hamiltonian to describe the single spectrum and wave functions. Due to the limited

Hilbert space of the tight—binding Hamiltonian, sum rule is severely violated (less

110

than one third of the dipole sum rule is accounted for). Since the RPA problem is
completely solved in the framework of tight—binding, one may be tempted to use this
method to study the dynamical responses of other clusters where a TB Hamiltonian

is available, or at least for a rough idea of the system under study.

There are also many ways the work presented in chapter 4 can be extended.
The most obvious one is probably to treat the electrons explicitly instead of being
treated implicitly through the use of a many—body—alloy (MBA) potential. There are
a couple of options one may take. The simplest and most computationally efficient
approach is to use a tight—binding type Hamiltonian for the valence electrons. The
most ambitious way is to use the ab initio methods like LDA. Although these ab initio
methods as of now are computationally not practical to use in simulations to extract
the density of states, which is the most useful statistical quantity as we did using
the MBA potential, they can be used in molecular dynamics simulations to study
other properties of finite systems. For example, one may use MD to generate a large
number of configurations for clusters at finite temperatures, and then to take these
configurations as input to study the optical properties like the plasmons. By using MD
to generate a canonical ensemble of configurations and taking the average of different
configurations, the effects of temperatures will be properly taken into account. For
example, the finite temperature will increase the width of the plasmon. The molecular
dynamics can also be used to study atom—cluster or cluster—cluster collisions. The
study of cluster collision may lead to new structural as well as dynamical properties.
It may also reveal how the the energy is transferred from the projectile to the cluster
target and how the energy is partitioned among the atoms in the target and how fast

equilibration is achived.

Appendix A

Multipolar Responses for C60

The particle-hole state from chapter 3 is reproduced here,

< I‘ll?h >= [£2 [fpri)fh(Ri)¢(2)(7‘i) + gp(Ril'g hR( )¢1( )] Y00(r )+

47, L2] WMIL ),gh _m( R) + fh(Ri)gp.——(mRn)7V)l¢O( (7‘1)1/1m(r1)+
1_07r _;(— mlgpr )X gh( Rz')lg,_.,,, ¢i(7‘i)Y2m(f‘i)
_ _Z Aph(z(A-1)
The quantities Aph ' )s are defined as
A2110) = éfpm) R:’)(th HYOO (Pi)
A2,,(i) = \/;EP(R') gthi )Yoofrz‘)
Aim = —Z(—4, >112.me 0 + f1(R.-)sp._m(R.-)1Y1m(h)
A2110) = “EX—1071’)® ghf Ri)lg,_m Y2m(f‘,-) (A2)

Note that there are two A2,,(2 ) s for l — 0, they correspond to the two (1)0(r i)’ s, that

is (1)0(r,) = ¢3(r,) and (1)003) = ¢1(r,) respectively. <I>l(r,-) and <I>2(r,-) are obviously

111

112

given by (1)1(7‘5) = ¢0(r,)¢1(r,) and <I>2(r,-) = dim). Now we consider the expansion
of the Coulomb interaction. First we represent the Coulomb interaction through its

Fourier representation

62 62

V —— / exp[ik- (r1 — r2 + R1 — R2)]dkdk (A.3)

 

: [Pl—P2+R1—R2] Z2772

Now we expand the plane waves into spherical harmonics using the following formula
expuk - r) = 4w 2 i’iz(kr)fim(f)ifl:.(l¥)- 01.4)

lm
After we carry out the integration of the product of four spherical harmonics, we

obtain the following formula for the Coulomb interaction

JM (L JM
I’L’

V Z Z {Z VIZJ’L’(7‘117'2)[Y1(IA‘1)® YL(R1)]JM [l/l'ff‘z) ‘3 YL'(R2)]* } (A-5)

with Wi,1/LI(T1,T2) defined by

 

 

32
v;i,,,L,(r1,r2) = 2.111W2H1X2L + 1)(21/ + 1)(2L’ + 1)
CiiisCflé’soz"+L"’-L’ [110mnuke.)ju(krs)js(kRs)dk (as)

We have used the definition for the tensor product

[Yidf‘ll ® ”203)le = Z Clifn/Illgmgyhmi(Elli/12m2(f.2)' (A7)

m1m2
The quantities inside the curl bracket define the multipolar expansion of the Coulomb
interaction. The quantum numbers J M define the multipolarity of the plasmons when
we calculate the response functions. Since C60 has spherical symmetry, all the matrix
elements are the same for a given J and different M ’s. The matrices using (1L) and

(l’ L’ ) as their indices for different J ’s are

113

J20 4x4
J=1 6x6
J22 7x7

Now we will find the RPA Green’s function 0. First G0 can be put in the following

form

2(5)!) — Eh) *12 .' _ L
(r, r; w) =2: Aph (i)(r(I>1, )(Ep _ 5h)2 _ (w + 2.77),) E A ph(j)<I>)2(rJ) (A8)

1”! 111' 121'

 

Now we proceed to solve the RPA equation for a particular multipole JM by expan-

sion.

 

Got/GO: ZZAjghu) )<1>,(

ph. [11

111'

J 1W

/ Z A*§.i.0><1>..<r.) Z Vijwsm, r.) [14.006 YL.(R >1

lgj L2L3

[13301) (8) YL,( (R), )]M Z ZAghln )<1>,, (rk)drjdrkdrjdrk

wP1h1 13k

2(5p1_

(5P1 _ 5’11)

 

6:)
01' 277? 22m P1h1( )14(T' )

 

2(5 —5h) 12L
—— —E E Apl,,i)(r<1>) ,p . E 13*1 2
Ph 111' 1( )(Ep - €132 - (90 +2702 1,1, M
13L3

 

2(5 -)e;,
D12L2J3L3 BlgL3 P1 1 A14 j)(p r A9
P1h1(5p1_€h1)2_(w+in)2)2p§l;l§ P1h1(-) 14(J) ( )

Byg— _ ZW/A )[14(1~.) <8 11031)]wa as, (A.10)

for fixed JM and D’ll’l’l2ng is defined by

DhLl’M’2 = f(D11(Ti1)(plg(ri2)Vlj]Ll,12L2(T117'2)d7‘1dr2 (All)

114

Note that some arguments of the summation notation 2 come before it, this is to

ensure each block in the formulas has a clear meaning. Next we consider

 

 

2(e —e),)
GOVGOVGoz A],,,i)(r<I>1 f” , . 19“?“
a???“ l( i)Efpr'éihlz-(WJHUV ph
Dl2L2,l3L3 Z Bl3L3 2(EP1 — Shl) *14L4 DlgLQ, ,13L3

h . .
P1h1 P1 1(5P1 — €h1)2 _ (w '1' ”ll2 plh 1

25 E}, *6 .
228125;: (“’2 2) . ,Ai....0)<1>i.<r.-)=

 

 

pghg 153' 5219—2-5hz) (W + “7)
Z Z Aplh i)((pl 2(813 — Sh) ' 8112th DlszJngg
Ph fli 1( “)5” _5h )2“ w + “7)2 p

 

2 e , - e , .
E13L3J4L4D14L4JsLs Z ZBlsLs ( p h l . Bp61h1(j )(I)16(T'j) (A12)

It
[hit] lej p1 1(6131 — €h1)2 — (90 +17])2

And E’lLlJQLQ is defined as

2(5 - 8h)
El1L1,l2L2 = Bllh 1L1 P . BleLQ A.13
a: (as - €02 - (w + 20)2 p” ( )

 

Now the structure of the expansion is clear, in the middle we have

1

 

 

 

 

Ueff=D—DED+DEDED—-~=Dl+ED (A.14)
So we have
2(6 "5h) 1 .
G: Ag}; 2)(D (7",) P . A* 2 _
Eh: 12,7“ ’1 (ep — 5,,)2 _. (w + z”)2 E ph(])
Z Z Ainlh z)(q)! (“)7‘ 2(5p — (Eh) 22 BlszU 12142 13L3
10" l11 1 2 (EP — ‘5th (w +l71)2 1212 Ueff
‘353
Z Z 8’31? 2““ “ 5’“) 14*" ,. on» (is) (A15)
mhl uj pl 1( 101—51102 - (w + in)2 P1 1 4

When we look the response to an external field F (r,), we need to calculate

h)=2/A]h( i)<I>1(r,) )F(r,~)dr,-df',~ (A.16)

115

Then the response is determined by

Fewer) )(rlF )>=):g(p,1)€p( if]; 5:)+,,,,,g*(p,h)—

 

 

2(5p—5h) 1L] lLlL
' B*11 U61 1122
;;§[g(ph) (Pp-ShV—(wi-WZ M U

12142

 

Z [B*:’21[’al;;( 2(5p1—8h1)*(P11h1):l (A17)

.9
Plh-l 5P1 ”- €h1)2—-—(w + “llz

We use I" to denote the real coordinate, R the ion coordinate, r coordinate relative to
the 10m, 1. e. r’ — —R + 1‘. Then the external field 1s denoted by r’ JY JM(r ’) (for .—I — 0,
we use r’2Y00(ri’)). In order to do the integration over 1', we rewrite the multipolar

field in terms of the ion coordinate R and the relative coordinate r.

 

 

I . . . 47rJ 2.] +1 _ ,
r JYJMfl") = RJYJMfR) +\/ ( 3 )RJ 1": ClinijJ—lll’f—ml

 

 

. . 2 J—1J2J—1 2J1
Y1m,(I‘)YJ_1,M_m1(R) + \/ 7f( ) ( 15 )( + )RJ—2r2

ZCZmlJ— 2114— m1Y2m1(f‘)YJ-2-M—m1(fi‘) + (A18)

The terms with l > 3 have no contributions to the response function since the maxi-
mum angular momentum of site transition is l = 2. For l = 0, F(r’) = r’2Y00(r’) and

the expansion is given by

(R2+r2+2R-r)=—1—R2+

.rI2YOO(rI) = m

1
\/4ir

1 2 2\/47r ~ .
r + Br Y (R)Y* (1') (A19)

 

 

Now the calculation of g(p, h) is relatively simple, and we will not write down the
formulas. In the Born approximation the differential cross section for inelastic electron

scattering is

(120 e2m 2 4p’
dads; = (‘57) P77; (91,01), (A20)

 

 

116

where p and p are the initial and final linear momenta of the electron, q = p — p,
is the transferred linear momentum, m the mass of the electron, w the energy of the
excited state and S(w, q) is the spectral function (or the dynamical structure factor)

of the cluster, which depends solely on the properties of the cluster and is given by
1
S(w,q) = —Im/exp(iq ~ r)G(r,r';w)exp(—z'q' r')drder'dQ'. (A21)
7r

In order to do the integration, we expand exp(iq- r) and exp(—iq- 1") into spherical

harmonics,
eprq-r) =4wzi’111qrmmm 3,161). (A221
1771.
exp1—2'q- r') = 4w Z1-i1’jz<qu')14:,(f'mm(<1> (A23)

1m

Now j,(qr)Y}m(i‘) acts as the external field in the calculation of the structure factor
as 'r’JY:11\/[(I:I) does in the calculation of the response function. Once again, we ex-

press j,(qr)Y1m(i‘) in terms of functions of coordinates of R and 7" using the following

 

 

formula,
. A 4 0° . m
JL(qr,)YLM(r) : W Z Zll‘Hz-iL\/(2ll +1)(212 +1)Clllzm1213ms
2L + 1 [112:0
j1,(qr)j1,(qR) [mm 61 191111)] W (11:24)

Now we can define a similar gL(p, h) for each multipole L, then S(w, q, L) is defined
by

2(5p — 5h)
(522 — Ehlz — (w + in)

 

5(w,q, L) = 29(1), h) 29*(10111) -
ph

 

229mm “5““) B*$hL‘U§}il’12L2

ph 111.1 (519 '_ 5h)2 "‘ (w + my
‘2L2
2 e — a .
2 3*1’1" ( p‘ h‘) 29091,121) (A.25)

 

mm mm (5P1 — €h1)2 _ (w + 277)

117

Since

L
2 4m‘LY5M1c‘1)4vr(—1)LYLM(é1) = 4w<2L+ 1) (A26)
Adz-L

we may omit 47riLYLM(€1) and the summation over m everywhere in the formulas, and

include a factor 47r(2L + 1) in our final formula for the cross section. The differential

cross section is given by

 

(120 62m 2 4' I Lm”
dad“) = (723—) 1% Z 47r(2L+1)S(w,q,L) (A27)
' L:O

We take Lmasr = 20 in our calculations.

 

 

Appendix B

Multipolar Responses for
C20,70,100

Once again the particle—hole state is reproduced here,

< PlPh >= J}: [fleilfh(Ri)¢filri) + gp(R1') ' gh(Ri)¢f(r.-)l YooIi‘i) +

1

im

1071' _Zm:(_ )X gh( Ri)]2,-—m ¢i(7‘i)Y2m(f‘i)
_ —2 Ag}: M10?) lm(ri)
(mi
The quantities Ag}: )’s are defined as

A22(R,~) = Jgfpfli R-)(th
A22(R,) = J;gp(Ri‘)) gh(R)

A;)'Z‘(R.) = (‘1m) ngprR R')9h,—m(Ri)+fh(Ri)gP1—m(Ri)l

AiflRi) = (-1)m -§;r-ngp( )Xgh( Ri)l2,_m

118

AGEI- 1'") lfP(R1)gh,-m( R1)+fthilgp.—m(Ri)l¢0("‘i)¢1("‘i)YIm(i‘i)+

(13.1)

 

 

 

119

Note that there are two Agg(R,-)’s for (lm) = (00), they correspond to the two ¢0(”1)’Sa
that is (1)002) = 45302) and <I>o(r,-) = gbflri) respectively. In order to keep the number
of indices in the formulas to a minimum, we did not introduce another index for this

distinction, (1)1(r) and (D203) are given by (1)1(7‘ )= ¢o(r,)¢1(r,-) and (D203) 2 4510“,).

Using the notation in eq. (B.1) for the particle—hole state < rlph >, we can put

the single particle Green’s function in the following form,

(,r r; w) =2 2 A2,? (7")YzmII'i)

ph "hm“

(\D

 

(5 —5) [*1m, 1* A
(. _;h)§_ :+m)2/l lm (Ri')¢z'(rw)Ypm/(r.v) (8.2)
up c

The matrix elements of the Coulomb interaction are computed through its Fourier

representation, but this time we treat R1 — R2 as a single vector instead of treating

R1 and R2 as two separate vectors.

62 62

v = lrl—r2+R1—R2|= 2 2/exp[ik (rl—r2+R1—R2)]dlsdf<

 

(4w) 310 2 ‘
2 -----—-—-2P/;llj1(k7'1))/1m(r1)ylm(k)ZZ—mjl’(kr2Yl’m’(r2)Yl'm'(k)

I’m’

:1 Jr(kR12)YLM(Rn)YE111(R)dkdk (13.3)
LM

After the integration over lc of the product of three spherical harmonics is carried

out, one obtains

V = Z BzmMmr ,Mm(i'1)Y,fm,(i'2)YLM(R12)
mlrlri’LM
f11(197‘1)jl'(kr2)jL(lez)dk, (13.4)

where Bffiflm, is defined as

 

, (21+ 1)(2L +1) ,m,
8117431,"; ,2 il+L- 132mg 2\J 47r(2l’+ 1) CmLoCllmLM.

120

Now we will consider the expansion of the integral equation G = Go — Go VG. First
we will consider

2(5). — 6).)
(8p - 5h)2 - (w + in)2

 

GOVGO— _ Z 2: Ag}: )(.14.,. (1.)

lmi
plhl limit!

Z Alillhmml R11 )ZBIIJX ,(12m2/(p11 T11)jll(kr11)¢12(ri2)j12(kri2)JL(kR112'2)

l1771111 LM
‘2'"2‘2

2(51)’ "' Eh’)
(51w-~ <'311I)(2R-)(w(+2'77)2

 

cm.dr..deLM(Ri.z-.)A;%7$’(R12)

 

A*S$I’(Rii)¢ll(7‘z) ')I ”Tn/(1‘17) Z; 2 Ag}: M))/lm( .)
p’h, l’lmI1I
2(5). — 5),) 1 1
A*ph1m1(1:{ R1 )Dlmi miA2mIz R1
(511-Shy-(<.a+277)2)2(1;;1 1 ’ 1112 22 ph( 2)
‘2m2'2
2 8 I — 5 I a: 'm’ * A
( P h ) A gr)“(Ri/)‘I>zi(ri))Yl,m,(rp), (B5)

 

(5p, - 5h!)2 - (w +217)2

where the quantities Dzlm,.1,)2m2.-2 are defined as follows

_ LM ‘ . .
Dzlmnmmig - Z Bz.m.,z.m.YLM(Rm2)
LM

/ ¢l1(ri1)j11(kri1)(D12(Ti2)j12(kri2)jLUCRiiiz)driidrizdk

Next we consider

 

A 2(5 - Eh)
_ Z 2 1m . . . 1"
GOVGOVGO -- h lmi Aph(Rz)¢l(rl)Yfm(rz)(€p _ 5h)2 _ (w + l")2

p
p’h’ l’m’i’

 

2(5 - 5), )
, A*11m1(R_ )Dl _ ’1 _ AIQTZQ(R_2 ) P1 1 .
LEE; ph 11 1mm 2m2221:4:41 pl 1 (5131 __ 5h1)2 _ (w + “7)2

 

2(5 ' — 512') I’m
A*’3m3(R.- )D. m . ,1... . A",m;‘(R.- ) P , A* 3:1.(11 )
isggs M1 3 3 33 i H M 4 (Er-511')?—(W+277)2
4m4'4

121

 

14,-114, =1;:A::::< 14.1141; ESTES)...

p’h’ l’m'i'

*llhmi
Z A (R41) 2 Dllmlil112m2i2E12m2'i2113m3i3D13m3i3114m4i4

(lmltl (277321?
‘4'"4'4 l:sms‘s
2(6p’ - Ehl)

AétZHRu ) 2A*iwu(Ri')‘1’z'(7‘i')YszI(i‘wl, (3-6)

 

(5p, - 51102 - (w +271)

where E12m2.2,13m3.-3 are defined as

 

2(8 —- 5}.)
l m P *1 m3
Ebmm [Maia _ 2 AM? 2( (Ri2)(5p — 602 — (w + 2'77le ,3. (R15).

If one continues to calculate the successive terms in the expansion of GOVG, one can
see that the expansion involves higher and higher powers of the matrix product of E

and D. Once the rule of the expansion is clear, G = Go — Go VG can be solved as

 

A 2(510 ‘ 5h)
G E. 1;. A13“ (4114,. (“)(sp—th—(wwn)?

 

 

 

(Iml. I
A*il,:"’<R 161101114 —Z 2 Air}: R-W (Mme)
Plh’ lllrnlu'l
2(5p _ (5h) 1)( 1 )
A*m1ml(R; D
(510 _ 5h)?“ (0.) + Z77)2)211§1:.1 1 + ED 11m1i1112m2i2
lzmz'z

m 2 5 ' _ 5 ' * ’m' * A

Aéih12(R.2 ) ( P h) A LIhI(Ril)¢ll('ril)Yme(ril). (B.7)

(Ep' - 5w)2 *- (w +271)2
From here on, the procedure to calculate the response functions is exactly the same as
that detailed in Appendix A for C60. Since we have to take into account the different
orientations of the nonspherical clusters, we outline how to average over the directions

of q in the following. 5(w, q) is given by

S(w,q) = $Im/exp(iq - r)G(r, r';w)exp(—iq- r')drdfldr'dfl’. (B.8)

 

 

122

This amounts to the following calculation

/Q1(7‘i)Y1m(f‘i) exp(iq ° 1")de0 = exp(iq ' R,)/(D1(Ti)Yim(f‘i)

exp(z'q- r;)drd§2 = exp(iq - Rf) / (P1(7"i)Y1m(f‘i)47T Z illjldqril

l’m’
l/zimr(f‘i)l/I'm'(€l)drd9 = 47Till/zm(€l)€XP(iq ' Ri)/<I>z('ri)jz(qr,-) (39)

It is obvious that only Yzm(Q)exp(iq - Bi) and Yfmlfi) exp(—iq- Rj) (we have two
plane waves in the formula for S(w,q)) depend on the directions of q. To average

over 61 is fairly simple and amounts to do the following integral,

/ Yzmfii) exp(z‘q- Ri)1/,rm,(<a) exp(—iq- Rj)dq =

/Y1m(él) [rm/(fl) Z47r'iLjL(qRij)YLMfill’Ev/(RUMQ =

L 1W

 

Z47riLjL<qRinYzM<Rt> CisioCm’M (8.10)

LIV!

 

(21+ 1)(2L +1)
47r(21’ +1)

Note that we have treated R0 = Rl- — RJ- as one vector.

"Illllllllllllllllllls