Studies of Nanoclusters’ Collision
Depositions on a Crystalline Surface and Graphene, and Negative Restitution Coeﬃcient of
Nanoclusters

Kuniyasu Saitoh
December 9, 2010

c Copyright by
Kuniyasu Saitoh
20111Acknowledgments
I would like to acknowledge Prof. Hayakawa who is my supervisor of PhD course of science in Kyoto
University. I am very grateful that he taught me many things about physics, research, and scientificmind. I really appreciate all the advice he gave me for writing this thesis. I also gratefully thank Prof.
Wada and all members of Prof. Hayakawa’s group for their fruitful advice and comments on research.
I also thank all staﬀs of Yukawa institute for theoretical physics (YITP) and Division of Physics andAstronomy, Graduate School of Science, Kyoto University for their kindness and hospitality.
I am grateful to Prof. Luding, Prof. Kuninaka, Prof. Otsuki, Prof. Hatano, and people in the
granular material community for their help and fruitful discussion. I also would like to acknowledge
Prof. Brilliantov and Dr. Bodrova for their collaboration on the part of this thesis and agreement toinclude the results in this thesis. I also thank the members of department of mathematics in Universityof Leicester for their hospitality, where the part of this thesis have been carried out. I would like tothank Prof. Paquette for his critical reading of a couple of our published papers.
Parts of the numerical computation in this work were carried out using computers of the YITP.
My research was supported by (i) the Global COE program “The Next Generation of Physics, Spunfrom Universality and Emergence” from the Ministry of Education, Culture, Sports, Science andTechnology (MEXT) of Japan and (ii) the Research Fellowship of the Japan Society for the Promotionof Science for Young Scientists (JSPS) and the Grant-in-Aid of MEXT (Grant No. 21.1958). In
particular, I thank the bilateral international exchange program (BIEP) which enable me to stay at
University of Leicester and collaborate with Prof. Brilliantov and Dr. Bodrova.
I am grateful to the YITP program, the Global COE program from MEXT, and the Society of
Powder Technology, Japan for giving me the chance to organize the international workshop titled”Recent Progress in Physics of Dissipative Particles - From fine powders to macroscopic behaviors of
granular particles -” which was held from November 24 to 26, 2010 at YITP and is strongly related
to the part of results of this thesis.
Finally, I would like to thank my wife, Megumi, for her much support to my research and understanding of my life.

2

Curriculum Vitae
Education
April 2008 -PhD student, Division of Physics and Astronomy, Graduate School of Science,
Kyoto University

March 2006

M. A. (Science), Kyoto University

April 2004 - March 2006 Master course, Division of Physics and Astronomy, Graduate School of
Science, Kyoto University
March 2004

B. S. (Science), Kyoto University

April 2000 - March 2004

Bachelor course, Faculty of Science, Kyoto University

Fellowship
April 2009 - March 2011

Research Fellow of the Japan Society for the Promotion of Science (DC2)

November 2009 - December 2009
icester

Visiting research, Department of Mathematics, University of Le-

Work Experience
April 2006 - March 2008

Employee, DeNA Co., Ltd.

3

Publications
1. Kuniyasu Saitoh, Anna Bodrova, Hisao Hayakawa, and Nikolai V. Brilliantov,”Negative Normal Restitution Coeﬃcient Found in Simulation of Nanocluster Collisions”,
Phys. Rev. Lett. 105, 238001 (2010) featured as Editors’ Suggestion.
2. Kuniyasu Saitoh and Hisao Hayakawa,”Motion of a free-standing graphene sheet induced by the collision with an argon nanocluster:Analyses of the detection and heat-up of the graphene”,
Phys. Rev. B 81, 115447 (2010).
3. Kuniyasu Saitoh and Hisao Hayakawa,”Simulation of depositions of a Lennard-Jones cluster on a crystalline surface”,
Prog. Theor. Phys. 122, 1081 (2009).
4. Kuniyasu Saitoh and Hisao Hayakawa,”Rheology of a granular gas under a plane shear”,
Phys. Rev. E 75, 021302 (2007).

4

Abstract
In this thesis, we show the results of the molecular dynamics (MD) simulations of the collisions
between nanoscale objects such as (i) nanocluster depositions on a crystalline surface, (ii) nanocluster
depositions on a free-standing graphene, and (iii) oblique impacts of two nanoclusters. The aim of
this thesis is to analyze the mechanics of nanoscale objects and to explain the results of the molecular
dynamics simulations in terms of the continuum mechanics.
At first, depositions of amorphous Lennard-Jones nanocluster on a crystalline surface are numerically investigated. From the results of the MD simulation, we found that the deposited nanocluster
exhibits a transition from multilayered adsorption to monolayered adsorption at a critical incident
speed. Employing the energy conservation law, we estimate the incident speed at which evaporation
occurs during the impact. We demonstrate the estimated value well agrees with the results of the MD
simulation. We also explain that the scaled adsorption parameter η ≡ Nadh /Ncls , where Nadh and Ncls
are the number of atoms adsorbed on the substrate and the cluster size, respectively, is proportional to
the square of the incident speed. The boundary shape of the adsorbed nanocluster depends strongly on
the incident speed, and some unstable modes grow during the impact. We also analyze the wettability
between diﬀerent Lennard-Jones atoms, and temperature dependence of a deposited nanoclusters.
At second, the motion of a free-standing graphene induced by nanocluster impact is investigated.
The graphene is bended by the incident nanocluster and a transverse deflection wave is observed. We
find that the time evolution of the deflection is semi-quantitatively described by the elastic beam theory. The time evolution of the temperature profile of graphene is also measured, and we demonstratethat the analysis based on the least dissipation principle well explains the result in the early stage ofimpact.
At third, the oblique impacts of nanoclusters are studied by the MD simulation and theoretically.
In the MD simulations, we explore two models – Lennard-Jones clusters and particles with covalently
bonded atoms. In contrast to the case of macroscopic bodies, the standard definition of the normal
restitution coeﬃcient yields negative values for oblique collisions of nanoclusters. We explain thiseﬀect and propose a proper definition of the restitution coeﬃcient which is always positive. We
develop a theory of an oblique impact based on continuum model of particles. A good agreement
between the macroscopic theory and simulations leads to the conclusion that macroscopic concepts ofelasticity, bulk viscosity and surface tension remain valid for nanoparticles of a few hundred atoms.

5

Contents
1 Introduction
1.1 Nanoscience and nanotechnology . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Nanocluster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Argon nanocluster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4 Silicon nanocluster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.5 Nanocluster deposition on a surface . . . . . . . . . . . . . . . . . . . . . . . . . .
1.6 Collision of two nanoclusters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.7 Charge-up of nanocluster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.8 Mechanics of graphene and carbon nanotube . . . . . . . . . . . . . . . . . . . . . .
1.9 Quantum eﬀect and the method of empirical potential . . . . . . . . . . . . . . . . .
1.10 The aim of this thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.11 Organization of this thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8
8
8
9
9
9
9
10
11
11
12
12

2 Depositions of a Lennard-Jones nanocluster on a crystalline surface
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Molecular dynamics simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.3 Numerical analysis of temperature dependence . . . . . . . . . . . . . . . .
2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Time evolution of impacts . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.2 Incident speed dependences of the scaled cluster size and cluster adsorption
parameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.3 Transition from partial wetting to perfect wetting of the deposited cluster . .
2.3.4 Morphology of the final configuration of the adsorbed atoms in the cluster .
2.3.5 Wettability between diﬀerent Lennard-Jones atoms . . . . . . . . . . . . . .
2.3.6 Temperature dependence of a deposited Lennard-Jones nanocluster . . . . .
2.4 Discussion and conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13
14
15
15
16
17
17
17

3 Depositions of an argon nanocluster on a free-standing graphene sheet
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Molecular dynamics simulation of the impact . . . . . . . . . . . . . . . . . . . . .
3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Time evolution of the deflection . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 Analysis of the deflection . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.3 Analysis of the heat-up . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27
28
28
29
29
30
32
33
33

6

19
21
22
22
23
25

4 Negative normal restitution coeﬃcient for nanocluster collisions
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 MD simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Modified definition of e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4 Theory of an oblique impact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35
36
36
39
40
41

5 Summary

43

A Tersoﬀ potential
A.1 Explicit form of the Tersoﬀ potential . . . . . . . . . . . . . . . . . . . . . . . . . .
A.2 Forces between interacting three atoms . . . . . . . . . . . . . . . . . . . . . . . . .

44
44
45

B Brenner potential
B.1 Explicit form of the Brenner potential . . . . . . . . . . . . . . . . . . . . . . . . .
B.2 Forces between interacting three atoms . . . . . . . . . . . . . . . . . . . . . . . . .

48
48
49

C Spin eﬀect on the restitution coeﬃcient of nanocluster
C.1 Spin-less approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C.2 Angular displacement α . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C.3 Derivations of Eqs.(4.4), (4.5) and (4.13) . . . . . . . . . . . . . . . . . . . . . . . .
C.4 Spin eﬀect on the restitution coeﬃcient . . . . . . . . . . . . . . . . . . . . . . . .

51
51
52
53
54

D Elastic modulus and surface tension of a nanocluster
D.1 Elastic modulus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
D.1.1 Isotropic compression . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
D.2 Surface tension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
D.2.1 Model of interacting two spheres . . . . . . . . . . . . . . . . . . . . . . . .
D.2.2 Model of compressed two spheres . . . . . . . . . . . . . . . . . . . . . . .

55
55
55
56
56
58

7

Chapter 1

Introduction
1.1 Nanoscience and nanotechnology
In the last decade, science and technology of nanomaterials are explosively developed [1–3]. The purpose of nanoscience and nanotechnology is to control nanoscale objects in order to innovate nanoscale
devices. Because the nanoscale objects are intermediate between single atoms and macroscopic materials, the properties of nanoscale objects are often peculiar [4] and many researchers investigate such
peculiar aspects of nanoscale objects, i.e. nanoclusters, graphene sheets, carbon nanotubes, etc.

1.2 Nanocluster
Nanocluster is an aggregate containing 10 - 106 atoms which are bonded by the interatomic forces[5–7]. If the total charge of nanocluster is neutral, the nanocluster is called neutral nanocluster. For
example, argon nanocluster is a typical neutral nanocluster, in which argon atoms are bonded by thevan der Waals interactions [8, 9]. On the other hand, gold nanocluster, in which constituent gold
atoms are bonded by metallic bonds and free electrons play an important role, is one of the moststudied metallic nanocluster [10–13]. Because nanoclusters are categorized by the species of atomand the interatomic interaction [14–17], and each aspect of nanocluster is extremely diﬀerent, we
restrict ourselves to the study of neutral nanoclusters such as argon nanocluster and silicon nanocluster.
Therefore, the eﬀect of charge, which we describe in Sec.1.7, is ignored in this thesis.
The equilibrium states of neutral nanocluster are characterized by temperature and size [18–31].
In the case of extremely low temperature, we can see face centered cubic (FCC), decahedral, andicosahedral structures. On the other hand, if the temperature is relatively high, we can see ”liquid”
nanocluster in spite of size [22, 25, 32]. Because the number of atoms in nanocluster is 106 at most,the finite size eﬀect plays an important role in the physics of nanoclusters. For example, thermal
fluctuation of each atom can influence the structure of nanocluster, and negative thermal conductivity
and smooth jump of heat capacity are found in the melting of nanoclusters [32–36]. The number ofatoms on surface of nanocluster is equally matched for the number of atoms inside the nanocluster.
Therefore, nanocluster has an extremely high reactivity and the van der Waals force is most importantto the interaction between two neutral nanoclusters [5, 7].
In experiment, nanoclusters are made by cooling down a supersaturated vapor [6, 37–39]. Thereare essentially two ways to cool down the supersaturated vapor: (i) employing an adiabatic expansion
through a nozzle into a high vacuum region, or (ii) mixing with a cold rare gas flow. Nanoclusters canalso be made by hitting a material with a laser pulse or an ion beam.
The nanoclusters produced in laboratory are detected by diﬀraction of electron beam [6, 37–39].
In this method, the fast electrons are scattered from the nanocluster and the diﬀraction rings aroundthe primary position of the electron beam is recorded. The interpretation of the diﬀraction profiles is8not straightforward, however, the size, structure, and temperature of nanoclusters can be detected bythis method.

1.3 Argon nanocluster
Argon nanocluster is one of the most studied nanoclusters and has been investigated experimentallyand numerically since the early 1980s [8, 9]. Because argon is a rare gas, it is a typical neutral
nanocluster. The equilibrium state of argon nanocluster depends on the size and temperature, andthe FCC, decahedral, and icosahedral structures are possible in the low temperature regime [22, 25,32]. However, the glassy state of argon nanocluster, i.e. amorphous argon nanocluster, was foundin the electron diﬀraction experiment [40], and this result triggered many numerical simulations ofamorphous Lennard-Jones nanocluster [41–45]. From the results of this experiment, the amorphousargon nanocluster is stable if the number of argon atoms is up to N = 800, however, if the number ofargon atoms becomes larger, the structure is changed to the multilayered icosahedron [40].

1.4 Silicon nanocluster
Silicon nanocluster is one of the most important nanocluster for practical purpose because of theintense photoluminescence at room temperature [46–49]. There are many experiments to produce
neutral silicon nanoclusters by hitting a laser pulse to a bulk silicon. Silicon nanoclusters can also
be produced by injecting vapor-phase silicon into an hydrogenated thermal plasma [49]. In this case,the surfaces of silicon nanoclusters, which are deposited on a substrate, are hydrogenated, and it
is expected that the surfaces of silicon nanoclusters are coated by hydrogen atoms during the free
flight. From the results of experiments and numerical simulations, it is known that silicon nanoclusterchanges its shape from prolate structure to spherical diamond structure around the number of constituent atoms N = 20 [46].

1.5 Nanocluster deposition on a surface
It is highly important to investigate the behavior of nanocluster depositions on a substrate for fabrication of high-quality films used in nanoscale electronic devices and photonic devices [50]. The ionizedcluster beam (ICB) technique makes it possible to control the translational kinetic energy of nanocluster [51–53]. Therefore, there have been many experimental and theoretical studies of nanocluster
depositions on a substrate [54–61].
The outcome of nanocluster deposition is mainly influenced by the incident speed of nanocluster[54, 55]. If the incident speed is extremely high, the cratering or implantation of the substrate isobserved [62–68]. On the other hand, if the incident speed is low enough, nanocluster is adsorbed
or rebounded on the substrate [69–72]. In this thesis, the incident speed is relatively low (0.01eV atmost), and the adsorption state of deposited nanocluster is mainly analyzed.

1.6 Collision of two nanoclusters
Collisions of two spherical particles are one of the most fundamental problem in physics, and it isalso important in chemistry, engineering, and many other fields [73–77]. The pioneering paper ”On
the contact of elastic solids” in 1882 by Hertz [78] analyzed contact mechanics of elastic bodies andproposed the Hertzian distribution of contact pressure which acts on the contact area between colliding
two macroscopic spheres. Because the origin of the Hertzian contact pressure is only elasticity, contact
forces caused by the other eﬀects, for example, the surface energy and viscoelasticity of macroscopic

9spheres, are investigated by many researchers [79, 80]. Especially, the JKR theory for the adhesive
forces is well known [81], and it is applied to many fields of granular physics, astrophysics, etc.
Although there are much eﬀort to understand the normal head-on collisions, the role of the tangential forces parallel to the contact plane during oblique collisions are not well understood. Indeed,
the process is aﬀected by the stick and slip motion, the rolling frictions and the eﬀect of roughnessof particle’s surface, and the mechanism of the oblique impact has not been well explained theoretically [79, 80].
Inelastic collisions, where a part of mechanical energy of colliding bodies transforms into heat,
are common in nature and industry. Avalanches, rapid granular flows of sand, powders or cereals
may be mentioned as pertinent examples [82, 83]. Moreover, inelastic collisions play an importantrole in astrophysical objects, like planetary rings, dust clouds, etc. An important characteristic ofsuch collisions is the normal restitution coeﬃcient e. The concept of a restitution coeﬃcient, as abasic one of the classical mechanics, has been introduced long ago by Newton; it addresses an impact
of macroscopic bodies. According to a standard definition, it is equal to the ratio of the normalcomponent of the rebound speed, g 0 (prime states for the post-collision value), and the impact speed,gg0 · n
.
(1.1)
e=−
g·n
The unit inter-center vector n = r12 /|r12 | at the collision instant (r12 = r1 − r2 ) specifies the impact
geometry. Since particles bounce in the direction, opposite to that of the impact, e is positive, e > 0,and since the energy is lost in collisions, e is smaller than one, that is, 0 ≤ e ≤ 1. This is a common
statement in the majority of mechanical textbook, where it is also claimed that e is a material constant.
Recent experimental and theoretical studies show, however, that the concept of a restitution coeﬃcientis more complicated: First, it depends on an impact speed [84–87], second, it can exceed unity for a
special case of oblique collisions with elastoplastic plate [88–90], where the energy of normal motion
can increase at the expense of the energy of tangential motion [88–90].
The increasing interest to nanoparticles, inspired by its industrial significance, raises an importantquestion, to what extent the macroscopic concepts are applicable and whether they acquire new features at a nanoscale. The collisions of nanoclusters has been studied in detail numerically [72,91–95].
It was observed that the surface eﬀects, due to the direct inter-cluster van der Waals interactions, play a
crucial role: The majority of collisions of homogeneous clusters, built of the same atoms, lead to a fusion of particles [91, 92]; they do not fuse for high impact speeds, but disintegrate into pieces [91, 92].
This complicates the analysis of restitutive collisions, which may be more easily performed for particles with a reduced adhesion. Among possible examples of such particles are clusters of covalently
bonded atoms, especially when their surface is coated by atom of diﬀerent sort, as for H- passivatedSi nanospheres [72]. These particles can rebound from a substrate, keeping their form after an impact
unaltered [72]. The bouncing nanoclusters demonstrate a surprising eﬀect – the normal restitution
coeﬃcient can exceed unity even for strictly head-on collisions [93].

1.7 Charge-up of nanocluster
The eﬀect of charges in ionized nanocluster is also an important issue of experimental and theoreticalstudies of nanocluster. Recent experimental technique of the FLASH free electron laser (FEL) is able
to ionize nanocluster by photo-absorption of atoms in nanocluster [96]. Because the energy of theelectron beam is much higher than the ionization potential (21.564eV, 8.152eV, and 13.598eV for anargon atom, a silicon atom, and a hydrogen atom respectively), the nanocluster is charged up during
the irradiation process, and finally the Coulomb explosion occurs.
It is also possible to charge up nanoclusters by a collision. In this case, the kinetic energy ofincident nanoclusters must be higher than the ionization potential of atoms in nanocluster. In this10thesis, the incident speed of nanocluster is set on a low value (0.01eV at most), thus the eﬀect ofcharge-up is neglected. However, the influence of ionization on collision dynamics is highly important
as a future work.

1.8 Mechanics of graphene and carbon nanotube
Graphene is a two-dimensional atomic layer of carbon atoms packed into a honeycomb lattice. Because graphene can be wrapped up into fullerenes and rolled into carbon nanotubes, graphene is
the most fundamental nanoscale carbon material and such flexibility of graphene encourages many
researchers to investigate its mechanical properties [97–101]. Because graphene is a monolayer ofcarbon atoms, its thickness can not be defined precisely in the sense of the continuum mechanics.
Therefore, the bending modulus
dM∼ Yh3 ,
(1.2)
D=dκ
where M, κ, Y, and h are the bending moment, the bending curvature, Young’s modulus, and thethickness of graphene, respectively, can not be defined precisely, too. However, the bending modulus
D can be measured experimentally by atomic force microscope (AFM) indentation, and diﬀerent
physical origins of D have been sought in experiments and numerical simulations [102–107]. Thereare essentially two ways to define the thickness h of graphene: (i) h = 0.34nm, which is the distance
between two carbon layers in graphite, and (ii) h  0.34nm, which is an eﬀective thickness mainly
proposed by the numerical simulations. Although there is an ambiguity in the definition of h, grapheneis thought to be an extremely rigid plate with Y = 1 ∼ 3TPa.
Carbon nanotube, on the other hand, is a quasi-one dimensional tube of carbon atoms. Manyexperimental and numerical studies have clarified the super rigidity of carbon nanotube [108–111],
and it is a promising candidate for the next generation rigid fibers.

1.9 Quantum eﬀect and the method of empirical potential
There are several methods for the numerical study of nanocluster. Thus, the choice of an appropriate
one depends on the size of nanocluster, the time scale, and the situation.
In the case that nanocluster contains only a few atoms, for example the number of atoms in nanocluster is in the range of 2 ≤ N ≤ 10, the first principle calculation, i.e. ab initio calculation ofquantum mechanics or Hartree-Fock method, is used to take into account the quantum eﬀects [112].
In this method, all electron calculations are possible [113], however, this method becomes cumbersome for larger size nanoclusters and the long time calculations are almost impossible.
The method based on density-functional theory [114–117] and quantum Monte Carlo simulation[118,119] can be accurate calculations, if the correlation terms, which is treated under approximations,are adequately calculated. Although the accuracy of approximations depends on the systems, thesemethods can deal with one hundred atoms by the semi-quantum approach.
In the case that the number of atoms is lager and the time scale is longer, classical interaction
potentials, which are based on approximate quantum models, are necessary. These interaction potentials, i.e. empirical potentials, contain parameters fitted by experimental results or the first principle
calculations. For example, the Lennard-Jones potential, which is a popular model for rare gases, contains the short-range repulsive part arisen from the wave function overlap between interacting two
atoms, and the long-range attractive part arisen from the electric dipole polarization of two atoms (theLondon dispersion force).
In this thesis, because the number of atoms in nanocluster is ranged between a few hundreds anda few thousands, the classical molecular dynamics simulations with empirical potentials are adopted.
Because we also use the Tersoﬀ potential [120–125] and the Brenner potential [126–129] for the11calculations of carbon and silicon atoms, respectively, the details of them are described in Appendix
A and B.

1.10 The aim of this thesis
In this thesis, we show the results of the MD simulations of the collisions between nanoscale objects
such as (i) nanocluster depositions on a crystalline surface, (ii) nanocluster depositions on a freestanding graphene, and (iii) oblique impacts of two nanoclusters. The aim of this thesis is to analyze
the mechanics of nanoscale objects with the help of the continuum model to explain the results of themolecular dynamics simulations.

1.11 Organization of this thesis
This thesis is organized as follows: In Chapter 2, nanocluster depositions on a crystalline surface are
numerically investigated. To clarify the adsorption state of deposited nanocluster, the critical incident
speed, at which the deposited nanocluster forms a monolayer, is investigated. From the analysis basedon the energy conservation law, the critical incident speed is estimated. In addition, the boundary
shapes of adsorbed nanocluster are analyzed. We also analyze the influence of the interaction between
diﬀerent species of atoms and temperature dependence of deposited nanoclusters [130]. In Chapter
3, nanocluster depositions on a free-standing graphene sheet are investigated. The deflection waves
propagating in the graphene sheet are observed, and the dynamic motion of the graphene is analyzed.
Employing the linear theory for elastic plate, the time development of the deflection wave is reproduced in the early stage of the impact. An anisotropic temperature profile of the graphene is also
observed, and such a profile is quantitatively explained by the least dissipation principle for an irreversible process [131]. In Chapter 4, oblique impacts of two nanoclusters are investigated. Because of
the reorientation of the contact surface and the finite contact time, the conventional restitution coeﬃcient becomes negative when the incident angle is large enough. The modified form of the restitution
coeﬃcient is introduced, which is always positive, and the macroscopic model for dissipative particles
well reproduces the results of the MD simulations [132]. In Chapter 5, all the results and conclusions
are summarized.

12

Chapter 2

Depositions of a Lennard-Jones
nanocluster on a crystalline surface

132.1 Introduction
A nanocluster containing 10 − 10, 000 molecules exhibits intermediate properties between bulk materials and individual molecules. Recently, there has been growing interest in the physics of nanoclusters. [4–7, 32] In particular, it is important to investigate the depositions of nanoclusters on solid
surfaces for the construction of high-quality films used in nanoscale electronic devices and photonic
devices. [50]The ionized cluster beam (ICB) technique was developed by Yamada et al. [51, 52] The ICB
technique is used to produce atomic clusters by employing adiabatic expansion of condensed vapour
through a nozzle into a high-vacuum region. In the ICB technique, clusters are ionized by electron
impact and then accelerated toward a substrate. Because the ICB technique controls the translational
kinetic energy of the cluster, there have been many experimental and theoretical studies aimed atunderstanding the eﬀect of the incident velocities of the cluster on the outcome of impact. [59–61]
The outcome of such a cluster impact is largely aﬀected by the incident velocity, as seen fromthe phase diagrams in Refs. [54, 55] If the translational kinetic energy per atom becomes too large,
the cluster can damage the substrate, and the cluster can break into pieces after the impact. [62–68]
However, if the translational kinetic energy per atom is less than 100 eV, the cluster is adsorbed on
the surface or reflected by the surface. Awasthi et al. carried out molecular dynamics simulations forcollisions of Lennard-Jones clusters with weakly attractive surfaces. [70,71] They discovered that thecluster rebounds when the translational kinetic energy of the cluster is larger than the adhesion energy.
Moreover, they clarified that a transition from adhesion to rebound occurs at the critical translational
kinetic energy. Järvi et al. carried out molecular dynamics simulations of low-energy deposition ofindividual metal clusters on a (100) surface. [69] They revealed that the heat generated by the impact
partially or completely melts the deposited cluster. As a consequence, the atoms in the cluster are
rearranged and adjusted to the atomic structure of the substrate. They found the maximum size at
which single clusters align epitaxially on the substrate.
Recently, Kuninaka and Hayakawa have carried out molecular dynamics simulations of two identical colliding clusters and investigated impact phenomena of nanoclusters subject to thermal fluctuations. [93] They found super-rebound events in which the restitution coeﬃcient is larger than 1.
They confirmed the validity of the macroscopic quasi-static theory of cohesive collisions. [86] Thissuggests that the research on nanoclusters is relevant even for the study for fine powders whose diameters range from 100 nm to 1 µm. [73, 74] They also revealed the mechanism responsible for thesuper-rebound process, the normal rebound, and the merging.
Although early numerical studies assumed that the clusters are highly crystallised, we also need toknow the properties of amorphous clusters. Indeed, it is easy to form metastable amorphous clusters
in terms of the quench process from high-temperature liquids. [41–45, 133]
The main purpose of our paper is to understand the behavior of the deposited amorphous LennardJones clusters on the crystalline surface at zero temperature. Here, we report on our molecular dynamics simulation of the depositions with the small incident energies per atom, which lie in the so-calledsoft-landing regime (0 − 2 eV). In addition, we report on the wettability between diﬀerent LennardJones atoms.
From the analysis of the final configurations of the deposited clusters, we find the existence of amorphological phase transition from the hemi-spherical droplet to the monolayer film at the criticalincident speed. The multilayered adsorption state is described on the basis of the energy conservation
law. Furthermore, we find that there are some unstable modes of the boundary shape of the depositedcluster.
This paper is organized as follows. Section 2 consists of three sections. In §2.1, we explain themodel of our numerical simulation. We explain our setup of cluster depositions in §2.2. In §2.3, we
explain how to analyze temperature dependence of deposited clusters numerically. Section 3 consistsof five sections. In §3.1, we show some time evolutions of impact processes. In §3.2, we explain
14

Figure 2.1: A snapshot of our simulation of a nanocluster deposition. The incident cluster contains
300 atoms that are bounded by the Lennard-Jones potential. The substrate consists of a single layer(120 × 120) of atoms on a triangular lattice.
how the cluster size and cluster adsorption parameter depend on the incident speed after the impact.
In §3.3, we discuss the transition from partial wetting to perfect wetting of deposited clusters. In§3.4, we explain the morphological change of the final configuration of adsorbed atoms in clusters.
In §3.5, we discuss the wettability between diﬀerent Lennard-Jones atoms. In §3.6, we analyze thetemperature dependence of deposited Lennard-Jones nanoclusters on a crystalline surface. In §4, we
discuss our numerical results and summarize the conclusion.

2.2 Molecular dynamics simulation
2.2.1 Model
To investigate the nanocluster depositions on a substrate, we perform a molecular dynamics simulation. Figure 2.1 displays a snapshot of our numerical simulation. Because we are interested in neutral
nanoclusters and substrates, the electrostatic interaction between atoms is not considered. We assume
that the potential energy of the interaction between two atoms can be described using the LennardJones (LJ) potential:!
 σ !12
σαβ 6 
αβ
−
,
(2.1)
U(ri j ) = 4αβ
ri j
ri j
where subscripts α and β specify the species of LJ atoms and ri j is the distance between two atoms
labeled by i and j. Here, αβ and σαβ are, respectively, the strength of the interaction and diameter ofthe repulsive core between α and β atoms. If α , β, we adopt the cross parameters of the LJ potential
using the Lorentz-Berthelot rule asσαβ =(σα + σβ )
,
2

αβ =

√

α β .

(2.2)

We mainly investigate the case that the cluster and substrate are constructed using the same atoms, A.
Here, we borrow the values of LJ parameters and mass of a typical inert gas, argon. Therefore, AA ,σAA , and the mass of an A atom mA are 1.65 × 10−21 J, 3.405 Å, and 6.63 × 10−26 kg, respectively.
[134–136] On the other hand, to study the eﬀect of the interaction energy between diﬀerent atoms ofthe cluster and substrate, we introduce a C atom as the constituent of the substrate in §3.5. We use
the values of LJ parameters and mass of carbon to the C atom, in which CC , σCC , and the mass of C
atom mC are 3.86 × 10−22 J, 3.354 Å, and 1.99 × 10−26 kg, respectively. [134–136] We should notethat the interaction energy between the cluster and substrate, AC , is several times weaker than AA . In
the following, we omit the subscripts of the LJ parameters of the interaction between the A atoms andwe adopt simplified notations , σ. We also adopt m as the mass of an A atom. We usep
, σ, and m as
the units of energy, length, and mass, respectively. Thus, the unit time is given by τ = mσ2 /.

15

Figure 2.2: Illustration of an LJ cluster formation consisting of 300 A atoms. (a) Initial configurationof atoms in the gas phase at T = 1.0. (b) A liquid cluster obtained from the quench into T = 0.5.
(c) An amorphous cluster obtained from the quench into T = 0.01.
We use a single-layer surface, which involves 120 × 120 atoms on a triangular lattice as thesubstrate with the periodic boundary condition. [137, 138] We set the lattice constant to 21/6 σαα(α = A, C) as the equilibrium distance between atoms. To avoid the destruction of the substrate,each atom of the substrate is also tethered to its equilibrium position by an elastic spring. In actual impacts of nanoclusters on substrates, the energy induced by an impact is relaxed to the internal motionof the atoms of the bulk region of the substrate. To represent such energy relaxation, we simply introduce the viscous force proportional to its velocity. The introduction of the viscous force has anotheradvantage to reduce the unrealistic boundary eﬀects. Indeed, if we do not introduce the viscous force,the acoustic wave would be transmitted across the boundary. The atom of the substrate at ri satisfies
the equation of motion
mα

X d
dri
d2 ri
eq
=
−
U(ri j ) − k(ri − ri ) − λ ,
2
dri
dt
dt
j

(2.3)

P
eq
where j is a summation over the interacting pairs i and j, and ri is its equilibrium position, andmα (α = A, C) is the mass of an atom of the substrate. Because we consider the substrate as a rigid√
surface, we use the spring constant k = 1.0 × 103 /σ2 . We simply adopt the value λ = 1.0 m/σ forthe coeﬃcient of viscosity.
We adopt the velocity Verlet method forp numerical integration of the equation of motion for each
atom with the time step dt = 1.0 × 10−3 mσ2 /. To reduce computational costs, we introduce thecutoﬀ length σcut = 3.0σ to the LJ potential, and we adopt the periodic boundary conditions in thehorizontal xy directions and the free boundary condition in the vertical z direction. It should be notedthat the viscous force is evaluated as the value at the previous time step for the numerical integrationof Eq. (2.3).

2.2.2 Setup
We make an LJ cluster by the temperature quench [139] into the metastable phase of LJ fluid. [140]We prepare 32, 108, 255, 300, 500, and 862 atoms in a periodic box and equilibrate at the temperature
T = 1.0 with the number density 0.05σ−3 in the gas state (Fig. 2.2(a)). It should be noted that the unitof the temperature becomes , because we set the Boltzmann constant to be unity.
To equilibrate the√gas at a specific temperature, we perform the velocity scaling until τ = 2000σ m/ for the relaxation
to a local equilibrium state. We have confirmed the equilibration of the total energy in the initial
relaxation process, and we quench the gas to T = 0.5. After an equilibration, a weakly boundedliquid cluster is formed (Fig. 2.2(b)) and is quenched to T = 0.01 to make it rigid. This two-stepquenching is adopted to form one cluster from an initial gas state. Indeed, if we quench the system
into T = 0.01 directly, many small clusters appear. After this equilibration, we obtain an amorphouscluster (Fig. 2.2(c)). We place the amorphous cluster at 10σ above the substrate and give the cluster
the translational velocity Vimp to make it collide against the substrate. It should be noted that the
16

z

Vimp

y

x

4 layers
40 atoms

40 atoms

Figure 2.3: (Color online) A snapshot of our simulation of a deposition of a Lennard-Jones cluster.
The incident cluster contains 862 atoms which are bounded by the Lennard-Jones potential. The
substrate is a (111)-terminated fcc surface consists of 4 layers of 40 × 40 triangular lattice of atoms.
amorphous cluster is metastable to maintain its shape within our observation time. The incident angleof the √cluster to the substrate normal is zero. The incident speed of the cluster ranges from Vimp = 0.1
to 5.0 /m.

2.2.3 Numerical analysis of temperature dependence
To analyze temperature dependence of deposited nanoclusters, we use a (111)-terminated fcc crystalline surface for the substrate. Figure 2.3 shows the setup of numerical simulation. The (111)terminated fcc surface consists of 4 layers of 40 × 40 triangular lattice. We set the lattice constantto 21/6 σCC . To avoid the destruction of the substrate, each atom of the substrate is tethered to its
equilibrium position by an elastic spring. Thus the atom of the substrate at ri satisfies the equation of
motionX d
d 2 rieq
U(ri j ) − k(ri − ri ),(2.4)
mC 2 = −
dri
dt
j
P
eq
where j is a summation over the interacting pairs i and j, and ri is its equilibrium position. Because
we consider the substrate as a rigid surface, we use the spring constant k = 1.0 × 102 /σ2 . We adopt
the velocity scaling method every 0.5τ to equilibrate the substrate at a specific temperature. The
temperature of the substrate is ranged from T = 0.01 to 0.4. We adopt the periodic boundary
conditions in the horizontal xy directions and the free boundary condition in the vertical z direction.

2.3 Results
2.3.1 Time evolution of impacts
Figures 2.4 and 2.5 display the time evolutions of the impact of the LJ√ cluster of 300 atoms on the
crystalline surface. Figures√ 2.4(a)–(d) represent the case of Vimp = 2.0 /m, and Figs. 2.5(a)–(d) are
the case of Vimp = 4.0 /m.
The incident cluster moves toward the substrate with its translational speed Vimp (Figs. 2.4(a) and2.5(a)), and hits the substrate (Figs. 2.4(b) and 2.5(b)). After the hit, the cluster is only deformed to
be a hemi-sphere (Fig. 2.4(c)) for the small incident speed. If the incident speed is, however, larger
than a critical value, the deposited cluster is split into many pieces (Fig. 2.5(c)). After the impact, thedeposited cluster is adsorbed on the substrate and settles into the final configuration (Figs. 2.4(d) and2.5(d)).
We observe that the impact process
and final configuration depend strongly on the incident speed√
Vimp . In the case of Vimp < 1.7 /m, no atoms can escape from the cluster√ during the impact. By
contrast, some atoms evaporate during the impact process for Vimp ≥ 1.7 /m. If the incident speed
is relatively low, the final configuration is a hemi-sphere on the substrate, as in the case of a partial
wetting of a liquid droplet on a dry surface. The deformation becomes larger as the incident speed
17Figure 2.4: Time evolution √of an impact of an LJ cluster of 300 atoms on the substrate, where theincident speed is Vimp = 2.0 /m. See the text for details.

Figure 2.5: Time evolution √of an impact of an LJ cluster of 300 atoms on the substrate, where theincident speed is Vimp = 4.0 /m. See the text for details.
√
increases. Above Vimp = 3.3 /m, the deposited cluster is completely split into√fragments, and theabsorbed atoms on the substrate form a monolayer coverage. Above Vimp = 4.5 /m, the depositedcluster bursts into fragments, and the absorbed coverage is no longer characterized by one cluster.
At the moment of impact, the temperature of the deposited cluster increases because the initial
kinetic energy is transformed into internal motion. Then, the temperature decreases owing to the heat
conduction into the bulk region of the material through the contact area. [69, 141, 142] The configuration of the deposited cluster is changed into an energetically favorable position during cooling.
Furthermore, the atomic structure of the deposited cluster is adjusted to the substrate. Finally, the18configuration is frozen because of the loss of heat.

2.3.2 Incident speed dependences of the scaled cluster size and cluster adsorption parameter

Figure 2.6: (a) A plot of the scaled cluster size ξ and (b) a plot of the cluster adsorption parameter η
as functions of the incident speed for the cluster of 300 atoms.

Figure 2.7: (a) Plot of the scaled cluster size ξ and (b) plot of the scaled cluster adsorption parameterη as functions of the incident speed.
In our simulation, the main cluster is detected using the clustering algorithm. [143] Following
Allen and Tildesley, we adopt the critical atom separation rC = 1.6σ. After the cluster settles into thefinal configuration, we represent Ncls as the number of atoms in the cluster. With the number of atoms
in the cluster before the impact N, we introduce the scaled cluster size:
ξ≡

Ncls
.
N

(2.5)

If ξ = 1, no atoms can escape from the cluster after the impact. On the other hand, if ξ < 1, some
atoms evaporate during the cluster impact.
Let us define an absorbed atom in the cluster if an atom at r in the cluster satisfies the relation|r − r s | < rC , where r s is the position of its nearest neighbor constituent of the substrate. Using thenumber of these adsorbed atoms Nadh , we can introduce the cluster adsorption parameter:
η≡

Nadh
.
Ncls

(2.6)

If η < 1, the cluster is regarded as a multilayered adsorption. However, if η = 1, the deposited cluster
is perfectly spread on the substrate, and it is a monolayered adsorption.
19

Figures 2.6(a) and (b) plot the incident speed
dependences of ξ and η for the cluster√of 300 atoms.
√
We find that ξ equals 1 below Vimp = 1.7 /m, but it decreases above
Vimp = 1.7 /m. On the√other hand, η increases
with the incident velocity below Vimp = 3.3 /m, but it is satisfied to η ' 1√
above Vimp = 3.3 /m. Figure 2.7(a) plots several results on ξ for N = 255, 300, 500, and 862, while
Fig. 2.7(b) is η − η0 for N = 32, 108, 255, 300, 500, and 862, where η0 is η at Vimp = 0. It seems thatη − η0 is independent of the size of clusters, while ξ exhibits weak size dependence.
How can we understand the behaviors in Figs. 2.6 and 2.7? During the impact, the temperature in
the cluster increases because the kinetic energy is transformed into internal motion. [69] We assume
that the energy flux to the substrate Φbulk during the impact is small, and the temperature becomesmaximum T max when the speed of the center of mass of the cluster becomes zero. Thus, the energy
conservation law can be written as1
3
3
2
mNVimp+ NT 0 ' NT max + ∆S ,
2
2
2

(2.7)

where T 0 is the temperature of the cluster before the impact, and ∆S is the change in the surface
energy. With the introduction of the surface tension γ, the height of the deposited cluster h, thecontact radius of the deposited cluster R, and the ratio φ = h/R, ∆S is given by∆S = γ(2πφR2 − 4πR20 ),(2.8)

where R0 is the mean radius of the cluster before the impact. By introducing the mean area fractionof the contact area ρadh = Nadh /πR2 , ∆S can be rewritten as
!
2Nadh
2
(2.9)
− 4πR0 .
∆S = γρadh
From Eq. (2.7), T max satisfies
T max = T 0 +2∆Sm 2
Vimp −.
3
3N(2.10)

Because the binding energy per atom in the cluster is roughly equal to , the evaporation takes place
at T max ' . In our simulation, T 0 is much smaller than T max , and the last term on the right-hand sideof
√ Eq. (2.10) is negligible for large N. Thus, the evaporation is considered to take place near V√imp '3/m. In Figs. 2.6(a) and 2.7(a), the scaled cluster size becomes ξ < 1 above Vimp = 1.7 /m,which is consistent with the above estimation. For the clusters with 32 and 108 atoms, ξ decreases
faster than the other cases. In such cases, we cannot ignore the last term on the right-hand side ofEq. (2.10).
During the impact, an evaporated atom carries away the volume energy uV , which is the potential
energy per atom and the kinetic energy 32 T max from the cluster. We assume that the internal energy ofthe deposited cluster decreases because of the energy flux to the bulk of the substrate. Therefore, after
the cluster settles into the final configuration, the energy conservation law can be written as1
3
2
mNVimp+ NT 0 = ∆S + Φbulk + (1 − ξ)N Ē,22(2.11)

where Ē = 32 T max − uV is the energy carried away by an evaporated atom. Here, (1 − ξ)N represents
the number of evaporated atoms. If the incident kinetic energy is not large, the number of atoms in thecluster is approximately preserved during the impact. Therefore, it is reasonable that the scaled cluster
size satisfies ξ ' 1. If we assume φ ' 1, the energy conservation law (Eq. (2.11)) can be simplified asmρadh 2ρadh  4πγR20 Φbulk 3 η=V +−+ T 0  ,(2.12)

4γ imp2γ  NN
2
20

where we have used Eqs. (2.6) and (2.9).
We use ρadh = 0.91, because the adsorbed atoms should match the hexagonal lattice on the substrate. The mean radius R0 of a cluster consisting of N atoms satisfies R0 = r0 N 1/3 , where we use
r0 = 0.68σ as a fitting parameter. The solid line in Fig. 2.6(b) is the theoretical prediction (Eq. (2.12)),
where the surface tension γ ' 3.5/σ2 and the energy flux per atom ϕbulk = Φbulk /N ' 1.7 are thetwo other fitting parameters.
The second term on the right-hand side of Eq. (2.12) can be written as
!
ρadh
3
2
−1/3
η0 (N) =4πr0 γN− ϕbulk + T 0 .
(2.13)
2γ
2
It is interesting that η − η0 (N) is independent of the cluster size. Figure 2.7(b) plots our numericalresults η − η0 (N) for η < 1, which support the validity of the theoretical prediction (solid line).

2.3.3 Transition from partial wetting to perfect wetting of the deposited cluster

Figure 2.8: (a) Plots of the height h (filled circle) and contact radius R (open square) and (b) a plot of
the contact angle θ as functions of the incident speed for the LJ cluster of 300 atoms.
Let us introduce the radius of the equimolar dividing surface (Gibbs Surface) [143]Z ∞1dρ(r) 22
R =−r drρadh 0dr(2.14)

as the contact radius of a deposited cluster, where ρ(r) is the area fraction of the contact area withradial distance from the center of mass of the adsorbed atoms in the cluster. We also define the cluster
height h as zmax − z0 , where zmax is the maximum vertical position in the atoms in the cluster, and z0 is
the minimum vertical position. Assuming a meniscus shape to the deposited cluster, we geometricallycalculate the contact angle θ.
Figure 2.8 displays h, R, and θ for the deposited cluster consisting of 300 atoms as functions of theincident speed. We observe that the cluster height√ h decreases and the contact radius R increases as theincident speed increases. Above Vimp = 3.3 /m, the height and contact angle become zero, which
implies that the deposited cluster becomes a monolayer film and is perfectly wetting on the substrate.
In this regime, the monolayer film is spread further and its boundary is partially chipped. Therefore,
the contact radius decreases. The clusters consisting of 108, 255, 500, and 864 atoms are also perfectly
wetting on the substrate at critical incident velocities. However, it should be noted that the transition
from multilayer film to monolayer film is only the morphological change of the deposited cluster. In
the case of the cluster consisting of 32 atoms, the number of adsorbed atoms is too few to define thewetting parameters h, R, and θ.

21

Figure 2.9: Configurations of the adsorbed atoms in a deposited cluster of 300 atoms for√ each incident
speed. Vimp equals (a) 0.5, (b) 1.0, (c) 1.5, (d) 2.0, (e) 2.5, (f) 3.0, (g) 3.5, and (h) 4.0 /m.

2.3.4 Morphology of the final configuration of the adsorbed atoms in the cluster
The boundary shape of the contact area depends strongly on the incident speed. To investigate themorphology of the boundary shape, we define the radial distance of the boundary r = f (ψ). Here,r and ψ are the usual radial and azimuthal coordinates. We take the origin to the center of massof the adsorbed atoms in the cluster. Moreover, we define a dimensionless variable g(ψ) for theboundary [144] asf (ψ) − R.
(2.15)
g(ψ) =RP
We also use its Fourier representations g(ψ) = n gn einψ , with the integer n = 0, ±1, ±2, · · · .
Figures 2.9(a)–(h) are the variety of the final horizontal configurations of the adsorbed atoms in
the cluster consistingof 300 atoms. We find that there are three phases in the boundary shape. Below√
Vimp = 1.5 /m, the boundary shape is grainy (Figs. 2.9(a), (b)), and |gn |2 has some peaks at√ highermodes. Thismaybecausedbythesmall
numberofadsorbed
atoms.
In
the
case
of
1.5
/m <
√
Vimp < 3.0 /m, the deposited cluster is uniformly spread on the substrate (Figs. 2.9(c), (d), √(e), (f)).
Thus, peaks√of |gn |2 vanish and their boundaries can be fitted by circles. In the case of 3.0 /m <Vimp < 4.0 /m, the deposited cluster becomes a thin film or a monolayer film. In this regime, the2boundary is partially chipped (Figs.
√ 2.9(g), (h)), where |gn | has intense peaks at some modes in our
simulation. Above Vimp = 4.0 /m, the deposited cluster bursts into fragments and the number ofadsorbed atoms is too few to define R. Thus, we cannot define g(ψ) in this regime.
The thermal fluctuation of a circular geometry step is estimated as h|gn |2 i = T/2πβRn2 from theequipartition of energy among the gn . [144, 145] Here, β is the step edge stiﬀness. In our simulation,
the thermal fluctuation √is estimated as h|gn |2 i ∼ 0.1, while the |gn |2 has peaks ranging from 1000 to3000 above Vimp = 3.0 /m. Therefore, these intensive peaks reflect on the growth of some unstablemodes of gn during the spread of the deposited cluster on the substrate.

2.3.5 Wettability between diﬀerent Lennard-Jones atoms
To investigate the eﬀect of the interaction energy on the outcome of impact, we also perform themolecular dynamics simulation in which the LJ parameters of A and C are used for the atoms in thecluster and the atoms of the substrate, respectively. Henceforth, we call this situation A-C case. The
Lorentz-Berthelot rule (Eq. (2.2)) is adopted to calculate the cross LJ parameters AC and σAC . Other
simulation settings are the same as those described in §2.2.
Figure 2.10 shows the cluster adsorption parameter η and contact angle θ as functions of theincident speed for the LJ cluster of 300 atoms. We also plot the data in the case of the LJ parameters
22Figure 2.10: (a) Plots of the cluster adsorption parameter η and (b) plots of the contact angle θ asfunctions of the incident speed for both A-C case (open square) and A-A case (filled circle).
of the A-A case. We observe that the cluster adsorption parameter η remains low even if the impact
speed is high (Fig. 2.10(a)). We also stress that any A cluster does not become a monolayer film inwhich θ becomes zero during the deposition onto a C surface. This is because the wetted state of theargon cluster is unfavorable on the crystalline carbon surface, which results from AC < . Thus, it isclear that not only the incident velocity, but also the choice of composites is important to determine
the final configuration of the deposited cluster.

2.3.6 Temperature dependence of a deposited Lennard-Jones nanocluster
Critical temperature√
The cluster moves toward the substrate by the incident speed Vimp = 2.0 /m and collides with thesurface. After the impact, the cluster is deformed and adsorbed on the surface, and the temperature of
the deposited cluster is relaxed to the surface temperature. We observe the system for 4000τ after thetemperature of the deposited cluster is relaxed.
To investigate the stabilities of the deposited clusters, we introduce the evaporation rate which isdefined by1ρs= ,(2.16)
τev τν
where ρs is the adatom area fraction and 1/τν is a probability per unit time for adatom to evaporate
from the surface [146]. Because the cluster is deposited on the surface in vacuum, the evaporation rate
1/τev equals zero if the deposited cluster is stable. From the results of our simulation, the evaporation
rate 1/τev equals zero if the temperature is below 0.3. By contrast, 1/τev is finite and the system is
not in a steady state if the temperature is larger than 0.3. Therefore, the deposited cluster is stable
only below the critical temperature T C ≡ 0.3.
Adsorption parameter and morphology of the deposited cluster
Figure 2.11 displays the temperature dependence of the adsorption parameter η of the deposited cluster. Each value in Fig. 2.11 is averaged over the last 2000τ time steps. The adsorption parameter η
is almost independent of the temperature and satisfies η < 1 for T < T C . It can also be seen that η
decreases as the cluster size increases.
Morphology of the deposited cluster
Figure 2.12 exhibits the snapshots of the deposited clusters in the steady state, and Fig. 2.13 displays
the temperature dependence of the scaled value R3 /N. Each value in Fig. 2.13 is averaged over the

23

[ε]
Figure 2.11: Plots of the adsorption parameter η of the deposited cluster as the functions of the temperature for N = 108 ( open squares ), 255 ( filled squares ), 500 ( open circles ) and 862 ( filled circles), respectively.

Figure 2.12: (Color online) Snapshots of the deposited cluster at equilibrium. The number of atoms
in the cluster is 500. The temperature of the cluster equals (a) 0.01, (b) 0.1, (c) 0.2 and (d) 0.3,
respectively.

[ε]Figure 2.13: Plots of the scaled contact radius R3 /N as the function of the temperature for N = 108 (open squares ), 255 ( filled squares ), 500 ( open circles ) and 862 ( filled circles ), respectively.
24

[ε]

[ε]

Figure 2.14: Plots of (a) the binding energy Eb and (b) the orientational bond order parameter Q6 ofthe deposited cluster as the functions of the temperature for N = 108 ( open squares ), 255 ( filledsquares ), 500 ( open circles ) and 862 ( filled circles ), respectively.
last 2000τ time steps. It can be seen that the contact radius R is almost constant for each cluster size,
and it seems that the morphology of the cluster does not depend on the temperature if T < T C .
Structural properties of the deposited cluster
To investigate the structural properties of the deposited cluster, we measure the binding energy and theorientational bond order parameter. The binding energy Eb is defined by the potential energy per atomin the cluster. The orientational bond order parameter is defined by a connecting pair of neighboring
atoms in the cluster. This parameter is related to the spherical harmonics function
Qlm ≡ Ylm (θ(r), φ(r)) ,(2.17)

where θ(r) and φ(r) are the polar and azimuthal angles of the bond r with respect to a fixed reference
frame. We also average Qlm over all bonds in the cluster as
Qlm =

1 X
Qlm (r),Nb r(2.18)

where Nb is the total number of bonds. Finally, to make this parameter invariant with respect torotations of the reference frame, we introduce the second-order invariants as the orientational bondorder parameter
v
u
t
l
4π X
Ql ≡
|Q |2 .
(2.19)
2l + 1 m=−l lm
It is known that Ql characterizes the dominant structure in the cluster [147–149].
Figures 2.14 (a) and (b) display the temperature dependence of Eb and Q6 , respectively. In our
simulation, we find that Eb and Q6 are steady even if the temperature T > T C . Each value in Fig. 2.14
is averaged over the last 2000τ time steps. It can be seen that Eb and Q6 increases as the temperature
increases below T C . When the temperature becomes larger than T C , Q6 decreases as the temperature
increases because of the evaporation.

2.4 Discussion and conclusions
In this paper, the incident kinetic energy per atom in the cluster is less than 2 eV. In this case, thedamage of the substrate due to the impact of a cluster can be ignored and we considered a single-layer
substrate. However, the eﬀect of the interaction between the deposited cluster and the bulk of thesubstrate is important. In general, the adsorption state is strongly aﬀected by the surface temperature,25

[150] but the substrate was assumed to be at T = 0 before the cluster impact in our simulation.
Therefore, the aﬀect of the surface temperature is also important for future study. Moreover, we
performed simulation only for one deposition event of the cluster at each incident speed. Thus, weshould take an ensemble average of impact processes for future study.
In conclusion, we found that deposited LJ clusters consisting 32, 108, 255, 300, 500, and 862
atoms exhibit a transition from√ multilayered adsorption to monolayered adsorption at the critical incident speed, Vimp = 3.3 /m. From our simulation, we clarified that the deposited clusters are
perfectly wetting on the substrate above the critical incident speed. Employing the energy conservation law, we estimated the critical value of the incident speed√at which the evaporation begins tooccur during the impact. The estimated critical value, Vimp = 1.7 /m, exhibits good agreement with
our simulation results of ξ. Using the energy conservation law, we also found that the scaled clus2 . These results
ter adsorption parameter is independent of the cluster size and is proportional to Vimp
exhibit good agreement with our simulation results.
√ Fourier analysis of g(ψ) and√ We performed thefound that some modes become unstable for 3.0 /m < Vimp < 4.0 /m. We also performed themolecular dynamics simulation of the A-C case and we found that the A cluster does not become amonolayer film on the C surface. Thus, we concluded that not only the incident speed, but also thestrength of the interaction between the cluster and the substrate is important to form a monolayer film
on a substrate. The deposited cluster is stable if T < T C , but the evaporation occurs if T > T C . The
morphology of deposited cluster is almost independent of the temperature for T < T C . The binding
energy of the deposited cluster Eb increases as the temperature increases. The orientational bond orderparameter Q6 increases as the temperature increases for T < T C , but decreases for T > T C because of
the evaporation.

26

Chapter 3

Depositions of an argon nanocluster on afree-standing graphene sheet27

3.1 Introduction
Graphene is a two-dimensional ( 2D ) atomic layer of carbon atoms on a honeycomb lattice. Recent
remarkable experimental techniques have made it possible to observe the motion of a free-standing or
suspended graphene sheet [98, 100]. Because electrons in a graphene can travel sub-micrometer distances without scattering, the study of graphene is active to make nanoscale electronic devices [101].
Graphene can be wrapped up into fullerenes and rolled into carbon nanotubes, and thus it is the mostfundamental structure of nano-carbon materials [97]. Such flexibility of graphene encourages many
researchers to investigate its mechanical properties. A recent experiment has detected the mechanical vibrations of suspended graphene sheets activated by radio frequency voltages, and has observed
vibration eigenmodes which are not predicted by the elastic beam theory [99]. In contrast to the electrical activations of graphene, it is also possible to activate the motion of graphene by nanocluster
impact [151]. The nanocluster impact can generate high pressure in localized areas of graphene, and
it is an appropriate method to verify the elastic theory for the plate deflected by the concentrated force.
In addition, nanocluster impact is also important for manufacturing nanoscale electronic devices on asubstrate [50, 55, 69, 130, 142]. Therefore, it is necessary to understand the motion of the graphene
induced by a collision with nanocluster in order to verify the elastic theory and to aim to construct thenanoscale electronic devices on a graphene sheet. However, there are a few studies which investigate
the time evolution of the local deformation of the graphene deflected by the nanocluster impact. In
this paper, we perform the molecular dynamics ( MD ) simulation to investigate the time evolution of
the deformation of a free-standing graphene sheet deflected by a collision with an argon nanocluster.
We find that analytic solutions of the elastic plate well reproduce the results of our MD simulation.
We also analyze the time evolution of the temperature profile of the graphene sheet.
The organization of this paper is as follows. In Section 3.2, we introduce our numerical modelof the nanocluster impact on a graphene sheet. Section 3.3 consists of three subsections. In Section
3.3.1, we show the time evolution of the deflection of the graphene. In Section 3.3.2, we analyze thetime evolution of the deflection. In Section 3.3.3, we analyze the heat-up of the graphene after the
impact. We discuss our results in Section 3.4 and conclude in Section 3.5.

3.2 Molecular dynamics simulation of the impact

Figure 3.1: (Color online) A snapshot of impact of an argon cluster on a free-standing graphene sheet.
The incident cluster contains 500 argon atoms. The graphene sheet contains 16032 carbon atoms on ahoneycomb lattice.
To study the dynamical motion of the graphene induced by a collision with anh argon cluster, wei
perform the MD simulation. We adopt the Lennard-Jones ( LJ ) potential φ(u) = 4 (σ/u)12 − (σ/u)6
for the interaction between two argon atoms with the distance u between two argon atoms, where weuse the LJ parameters [152]  = 1.03 × 10−2 (eV) and σ = 0.340 (nm). We also adopt LJ potential
for the interaction between an argon atom and a carbon atom, where we use the cross parameters√of LJ potential int and σint , which are defined by the Lorentz-Berthelot rule as int =  0 and
σint = (σ + σ0 )/2, respectively. Here,  0 = 2.40 × 10−3 (eV) and σ0 = 0.335 (nm) are the LJ
28

parameters for carbon [153, 154]. Finally, we adopt the Brenner potential, which is widely used for
simulations of a graphene and a carbon nanotube, for the interaction between two carbon atoms [126].
Figure 3.1 displays a snapshot of our impact simulation. The graphene involves 16032 carbon
atoms on a honeycomb lattice. The bond length of the graphene is 0.142nm and the length of one edge
is approximately equal to 20nm. The carbon atoms on the edges parallel to the x-axis are arranged inarmchair geometries, and the carbon atoms on the edges parallel to the y-axis are arranged in zigzag
geometries [155, 156]. The boundary conditions of the four edges of the graphene are free, and theinitial temperature of the graphene is 1.2K. The cluster containing 500 argon atoms is made fromargon gas by the temperature quench method [133, 157]. At first, we prepare 500 argon atoms in aperiodic box and equilibrate at 119.6K with the number density 1.27nm−3 in the gas state. We quench
the temperature to 59.8K. After an equilibration, a liquid-like argon cluster is formed. We further
quench the temperature to 1.2K to make it rigid, and an amorphous argon cluster is formed [130].
The center of mass of the amorphous argon cluster is placed at 5.1nm above the center of mass of thegraphene. The argon cluster is translated with the incident velocity V to collide with the graphene.
The incident angle of the argon cluster to the graphene normal is zero.

3.3 Results
3.3.1 Time evolution of the deflection

Figure 3.2: The deflection of the graphene sheet ζ at (a) 2.2 ps and (b) 2.8 ps after the initial hitting.
The incident cluster contains 500 argon atoms, and the incident speed is 316m/s.

Figure 3.3: The deflection of the graphene sheet ζ at (a) 2.2 ps and (b) 2.8 ps. The incident cluster
contains 500 argon atoms, and the incident speed is 790m/s.
Let us demonstrate the motion of the graphene induced by the collision with the argon cluster inthe case of V = 316 (m/s). Figures 3.2 display the time evolution of the deflection of the grapheneζ as a function of x and y coordinates. In this figures, we divide the xy plane into 32 × 32 cells andaverage over z-components of the positions of carbon atoms in the center-of-mass frame. We definet = 0 as the time at which the argon cluster contacts the graphene sheet. At the impact, the circular29

region around the center of the graphene is bended by the incident argon cluster ( Figs. 3.2 (a) ), andthe transverse deflection wave is isotropically propagated in the graphene ( Figs. 3.2 (b) ). In thelaboratory system, the graphene is moved downward and immediately reaches the uniform motion
along the z-axis with the speed 28.4m/s. During the impact, the incident argon cluster adsorbs on
the graphene and does not rebound. Figures 3.3 display the time evolution of ζ for the incident speed
V = 790 (m/s). At the impact, the circular region around the center of the graphene is strongly bended
by the incident argon cluster ( Figs. 3.3 (a) ), and the transverse deflection wave is observed ( Figs.
3.3 (b) ). During the impact, the incident argon cluster bursts into fragments and some fragments
are scattered and the rest of fragments adsorb on the graphene. We have also examined the cases of
V = 158, 474, and 632 (m/s), and the bending formation and the propagation of transverse deflection
wave are also observed. In all cases, the deflection wave in the graphene passes through the boundary
without reflection, and the graphene ripples after the impact. We have never observed any defectformations in the graphene sheet.

3.3.2 Analysis of the deflection
Equation of motion
To analyze the result of our simulation, we examine the linear theory of the elasticity in description of
the deflection of the graphene [158, 159]. Because the elastic properties of a 2D hexagonal structure
are isotropic [158], we ignore the anisotropic properties of the graphene sheet. Thus, the equation of
motion for the deflection is given byρζ̈(r, t) +h3 E
∆2 ζ(r, t) = p(r, t) ,12(1 − µ2 )(3.1)

Here, ρ = 7.59 × 10−7 (kg/m2 ) is the mass per unit area of the graphene, and ζ̈(r, t) represents∂2 ζ(r, t)/∂t2 . Because graphene is a single atomic layer of carbon, its thickness h is sometimes set tobe the diameter of a carbon atom, 0.335nm. However, Yakobson et al. indicated that h = 0.066 (nm)should be used in their simulation of single-walled carbon nanotubes [108]. We still do not have any
consensus on the proper value of h [102–105, 110]. Thus, to avoid ambiguous definition of the thickness, we use the thickness and the elastic moduli which are directly obtained from the analysis of theBrenner potential. Following Ref. 26, we use the thickness, Young’s modulus, and Poisson’s ratio as
h = 0.0874 (nm), E = 2.69 (TPa), and µ = 0.412, respectively. The right hand side of Eq. (3.1) is theexternal pressure due to the argon cluster impact. Because the deflection is symmetric with respect tothe z-axis, we assume that ζ and p depend on time t and the distance from the z-axis r.
Hertzian contact pressure
Although the external pressure changes during the impact, we simply assume that p(r, t) is an impulsive pressure with the distribution of the Hertzian contact stress. Thus, we may assumer r 23Fp(r, t) = −1−
Θ(a − r)δ(t) ,
(3.2)
a
2πa2
where F and a are the impulse and the contact radius of the incident argon cluster [158]
"

(
)#1/3
3FR 1 − µ2 1 − µ02
a=
+
,
4
E
E0

(3.3)

respectively. Here, δ(t) is Dirac’s delta function, and Θ(a − r) is the Heaviside function which isdefined as Θ(a − r) = 1 for r < a and Θ(a − r) = 0 for r > a. The mean radius, Young’s modulus, andPoisson’s ratio of the argon cluster are R = 1.6 (nm), E 0 = 3.69 (GPa), and µ0 = 0.396, respectively,
30

Figure 3.4: (Color online) The MD simulation results of the mean deflection of the graphene ( open
circle ) which are averaged over the azimuthal coordinate, and the solutions of the equation of motion,i.e. Eqs. (3.6) ( red solid line ) and (3.9) ( green broken line ) at (a) 2.2 ps and (b) 2.8 ps. The incident
cluster contains 500 argon atoms, and the incident speed is 316m/s. The magnitude of the impulse is
1.96 × 10−10 N · ps.
which are estimated from our another MD simulation [160, 161]. In addition, we assume that thecontact area of radius a moves downward with the speed V at the impact. Thus, the initial conditionsof the deflection are ζ(r, 0) = 0 and ζ̇(r, 0) = −VΘ(a − r). Because we consider the behaviors in thevicinity of the center of the graphene, we solve Eq. (3.1) as if the graphene sheet is infinitely large.
The Fourier transform and the Laplace transform of Eq. (3.1) yield
H(k),ζ̂k (s) = − 2
s + Dk4(3.4)

where D ≡ h3 E/12ρ(1 − µ2 ). Here, we introduce the function
H(k) ≡

3F sin(ak) − ak cos(ak)J1 (ak)+ 2πaV
,
3
ρ
k
(ak)

(3.5)

where Jn (x) is the Bessel function
represent the Laplace transform and the
R ∞ for an integer n. Here,R we∞−stFourier transform as ζ̂k (s) ≡ 0 ζk (t)e dt and ζk (t) ≡ −∞ dkζ(r, t)e−ik·x , respectively. The inverse
Laplace transform and the inverse Fourier transform of Eq. (3.4) yield√Z ∞sin( Dk2 t)ζ(r, t) = −H(k)J0 (kr)
dk .
(3.6)
√
0
Dk
A flat punch pressure
If we adopt a flat punch impulsive pressurep(r, t) = −FΘ(a − r)δ(t)
πa2

(3.7)

instead of Eq. (3.2), the Fourier transform and the Laplace transform of Eq. (3.1) yieldζ̂k (s) = −2H0 J1 (ak)1,aks2 + Dk4

where we introduce the constant H0 ≡ ρ−1 F + πa2 V. In this case, the solution of Eq. (3.8) is
√
Z
2H0 ∞
sin( Dk2 t)ζ(r, t) = −J0 (kr)J1 (ak) √
dk .
a 0
Dk2

31

(3.8)

(3.9)

Figures 3.4 display the time evolution of the deflection ζ(r, t) in the case of V = 316 (m/s). In
this figures, the open circles are our MD simulation results which are averaged over the azimuthal
coordinate. The red solid and green broken lines represent Eqs. (3.6) and (3.9), respectively. Here, we
use F = 1.96 × 10−10 (N · ps) for both Eq. (3.6) and Eq. (3.9). We have also examined the deflectionof the graphene in the case of V = 158 (m/s), and we find that the time evolution of ζ(r, t) is welldescribed by Eqs. (3.6) and (3.9) with F = 1.25×10−10 (N · ps). However, in the cases of V = 474, 632
and 790 (m/s), Eqs. (3.6) and (3.9) are no longer applicable because the incident argon cluster bursts
into many fragments which collide with the graphene, and the distribution of the external pressure canneither be approximated by the Hertzian contact stress nor a flat punch pressure.

3.3.3 Analysis of the heat-up

8
6
4
2
0

8
6
4
2
0

Figure 3.5: The temperature profile of the graphene sheet T (x, y) at (a) 2.2 ps and (b) 2.8 ps after the
initial hitting. The incident cluster contains 500 argon atoms, and the incident speed is 316m/s.
To study heat up of the graphene, we introduce the local temperature. We divide the graphene into
64 × 64 cells along the x- and y-axes and define the temperature of the j-th cell as
2m X
vi − u j ,Tj =3kB N j i=1
Nj(3.10)

where kB and m = 1.99 × 10−26 (kg) are the Boltzmann constant and the mass of carbon atom,
respectively. In Eq. (3.10), N j , vi and u j are the number of carbon atoms in the j-th cell, the velocityof the i-th carbon atom which is in the j-th cell and the mean velocity of the j-th cell, respectively.
The mean velocity of the j-th cell is defined as
1 X
vi .
uj =
N j i=1
Nj(3.11)In order to take a sample average of T j , we rotate the nanocluster around the line before the impact,
where we use the diﬀerent angle for each sample. Here, the line is parallel to the z-axis and intersects
at the center of mass of the nanocluster. If we project the 64 × 64 cells to the xy plane, Eq. (3.10)
approximately represents the temperature profile T (x, y). Figures. 3.5 display the time evolution ofT (x, y) which is averaged over 20 samples in the case of V = 316 (m/s). Although the thermal
conductivity of a 2D hexagonal structure is isotropic [162], the results of T (x, y) are anisotropic.
Let us explain the anisotropic profile of T (x, y). The nanocluster collides with the graphene in thevicinity of the center of mass of the graphene (X, Y). By the impact, the vicinity of (X, Y) is heated
up and the heat current q flows from (X, Y) to the edge of the graphene. Then, q is symmetrical withrespect to (X, Y), and we adopt (X, Y) for the origin. In such an irreversible process, thanks to the least32

dissipation principle, the rate of the entropy productionZ
D = − κ−1 q2 dA(3.12)

A

is expected to be minimum, where A and κ are the area of the graphene and the heat conductivity
per unit area of the graphene, respectively [163]. If we assume that κ is a constant, the variation
δD = 0 leads ∇ · q = 0 [163]. Therefore, from Fourier’s law of heat conduction, the deviation of thetemperature δT = T (x, y) − T 0 satisfies Laplace’s equation ∆δT = 0. Here, T 0 is the temperature of thegraphene before the impact. Because δT is finite at (X, Y), the general solution of Laplace’s equationis∞X
δT (r, θ) = const. +rn an cos(nθ + φ)
(3.13)
n=1

in the polar coordinate, where an and φ are the integral constants [164]. Because q is symmetricalwith respect to (X, Y), the integer n satisfies cos (n(θ + π) + φ) = cos(nθ + φ) and sin (n(θ + π) + φ) =sin(nθ + φ). Thus, n should be even. Therefore, δT (r, θ) is distributed around (X, Y) asδT (r, θ) = const. +∞X

r2m a2m cos(2mθ + φ) .

(3.14)

m=1

In Fig. 3.5 (a), the heated region can be seen as a quadrupole distribution around (X, Y) which isthe case of m = 1 in Eq. (3.14). On the other hand, in Fig. 3.5 (b), the heated region is no longer
distributed as Eq. (3.14). In this case, it seems that the least dissipation principle is no longer correct,
and it is necessary to solve the heat equation with the boundary conditions correctly.

3.4 Discussion
Although Eqs. (3.6) and (3.9) seem to well describe the results of our MD simulation, the solution does
not satisfy the boundary conditions, and these are not applicable except for the case that deformations
are localized in the vicinity of the center of the graphene sheet, i.e. immediately after the impact.
Note that it is diﬃcult to obtain an analytic solution of Eq. (3.1) which satisfies the completely freeboundary conditions [165, 166]. If we simply estimate the magnitude of the impulse from the change
in momentum of the incident argon cluster, F = 3.59 × 10−10 (N · ps) which is about two times
larger than the fitted value in Figs. 3.4. However, the value is over-estimated, because the loading
force can change during the impact and the dissipative force plays important role for the collision ofclusters [87, 93]. In the case of V ≥ 400 (m/s), the impact processes are further complicated by many
fragments of the argon cluster, and Eqs. (3.6) and (3.9) are no longer correct. Therefore, it is necessary
to improve the functional form of p(r, t). Note that if we use h = 0.335 (nm) in Eqs. (3.6) and (3.9),
the wave propagates much faster than the actual propagation observed in our MD simulation. Thus,
the thinner thickness h = 0.0874 (nm) is more appropriate. The analysis based on the least dissipation
principle reproduces our simulation result of the temperature profile in the early stage of impact.
However, in order to describe the time evolution of the temperature profile, it is necessary to solve theheat equation with appropriate boundary conditions.

3.5 Conclusion
In conclusion, we perform the molecular dynamics simulation of the graphene sheet induced by acollision with an argon nanocluster, and the bending formation and the propagation of transverse
deflection wave are observed. We find that the linear theory of the elasticity well explains the time
evolution of the deflection of the graphene, where the deflection is represented by using the analytic
33

expressions Eqs. (3.6) and (3.9). In addition, we conclude from the analysis of the motion of thegraphene that the actual thickness is much thinner than the diameter of a carbon atom. We also
analyze the time evolution of the temperature profile, and find that the analysis based on the leastdissipation principle reproduces our simulation result in the early stage of impact. We believe that thepredictions of the bending formation and propagation of transverse deflection wave are necessary forthe construction of the nanoscale electronic devices on a graphene sheet.

34

Chapter 4

Negative normal restitution coeﬃcient fornanocluster collisions35

4.1 Introduction
Inelastic collisions, where a part of mechanical energy of colliding bodies transforms into heat, arecommon in nature and industry. Avalanches, rapid granular flows of sand, powders or cereals may
be mentioned as pertinent examples [82, 83]. Moreover, inelastic collisions define basic propertiesof astrophysical objects, like planetary rings, dust clouds, etc. An important characteristic of suchcollisions is the so-called normal restitution coeﬃcient e. According to a standard definition, it isequal to the ratio of the normal component of the rebound speed, g 0 (prime states for the post-collision
value), and the impact speed, g
g0 · n
e=−
.
(4.1)
g·n
The unit inter-center vector n = r12 /|r12 | at the collision instant (r12 = r1 − r2 ) specifies the impact
geometry. Since particles bounce in the direction, opposite to that of the impact, e is positive, e > 0,and since the energy is lost in collisions, e is smaller than one, that is, 0 ≤ e ≤ 1. This is a common
statement in the majority of mechanical textbook, where it is also claimed that e is a material constant.
Recent experimental and theoretical studies show, however, that the concept of a restitution coeﬃcientis more complicated: First, it depends on an impact speed [84–87], second, it can exceed unity for a
special case of oblique collisions with elastoplastic plate [88–90], where the energy of normal motion
can increase at the expense of the energy of tangential motion [88–90]. Still, it is believed that e ≤ 1
for a true head-on collision.
The concept of a restitution coeﬃcient, as a basic one of the classical mechanics, has been introduced long ago by Newton; it addresses an impact of macroscopic bodies. The increasing interest to
nanoparticles, inspired by its industrial significance, raises an important question, to what extent themacroscopic concepts are applicable and whether they acquire new features at a nanoscale. The collisions of nanoclusters has been studied in detail numerically [72, 91–95]. It was observed that the surface eﬀects, due to the direct inter-cluster van der Waals interactions, play a crucial role: The majority
of collisions of homogeneous clusters, built of the same atoms, lead to a fusion of particles [91, 92];
they do not fuse for high impact speeds, but disintegrate into pieces [91, 92]. This complicates theanalysis of restitutive collisions, which may be more easily performed for particles with a reduced
adhesion. Among possible examples of such particles are clusters of covalently bonded atoms, especially when their surface is coated by atom of diﬀerent sort, as for H- passivated Si nanospheres [72].
These particles can rebound from a substrate, keeping their form after an impact unaltered [72]. The
bouncing nanoclusters demonstrate a surprising eﬀect – the normal restitution coeﬃcient can exceed
unity even for strictly head-on collisions [93].
In this Letter we investigate the oblique impact of nanoclusters with the reduces adhesion bymeans of Molecular Dynamics (MD) and theoretically, using concepts of continuum mechanics. Unexpectedly, we have found that the normal restitution coeﬃcient, as defined by Eq. (4.1), acquires forlarge incident angles negative values, e < 0. We explain this eﬀect by the reorientation of the contact plane during an impact and quantify it. Moreover, we propose a modified definition of e, which
preserves its initial physical meaning and yields always positive values. To describe the collision of
nanoclusters we develop a continuum theory. Surprisingly, the macroscopic approach quantitatively
agrees with MD even at nanoscale.

4.2 MD simulations
We study two models - a simplified model (A), which mimics interactions of nanoclusters with the
reduced adhesion and realistic model (B) for interaction of nanoclusters with covalently bonded atoms– H-passivated Si nanospheres. For the modelh A, which is less computationallyexpensive, we adopt
ithe Lennard-Jones ( LJ ) potential φ(r) = 4 (σLJ /r)12 − (σLJ /r)6 for the interaction between two

36

Figure 4.1: (Color online) Initial (left) and final stage (right) of the nanocluster collision. The initial
relative velocity is v12 (0) = V and the incident angle is γ. The unit normal n specifies the orientationof the contact plane. For large γ a noticeable reorientation of this plane is observed. Here the collisionof H-passivated Si nanospheres (model B) is shown.
hiatoms in the same cluster and the modified LJ potential φint (r) = 4 (σLJ /r)12 − c(σLJ /r)6 for theinteraction between atoms in two diﬀerent clusters. Here the cohesive parameter c = 0.2 controls theadhesive force [93–95] between clusters, while , σLJ , and r are, respectively, the depth of potential
well, the diameter of the repulsive core, and the distance between two atoms. The nanoclusters of
N = 500 atoms were prepared by the two step temperature quench to obtain a rigid amorphous
particle [130]. The diameter of nanocluster d was defined as the maximum distance between thecenter of mass of the nanocluster and the atom on the surface, so we find d = 10.46 σLJ . For themodel B we adopt the Tersoﬀ potential [120–125] for the covalent Si-Si, Si-H, and H-H bondings.
The Si nanospheres, containing 2905 Si atoms arranged in a diamond structure, are fully coated by
852 H atoms. The radius of Si nanosphere is about d = 2.4 nm.
We fix the√ modulus of the relative inter-cluster velocity v1 (0) − v2 (0) = v12 (0) = V and set it
to V = 1.0 /m and 1850 m/s for the model A and B, respectively. We vary the incident angle γ
between n and V (see Fig. 4.1), so that the normal impact velocity, Vn = V cos γ is changed. The
nanoclusters do not rotate before an impact and have zero angular velocities, ω1 (0) = ω2 (0) = 0. To
make an ensemble average, we randomly turn one of the clusters around the axis, passing throughits center and perpendicular to the contact plane. Due to rough atomic surfaces of the clusters, thisresults in varying contact configurations at each impact. Hence, for every incident angle γ we performaveraging over 100 collisions with diﬀerent contact conditions for model A and over 10 collisions
for model B. The clusters’ deformation during an impact is quantified by the normal displacement,ξn (t) = d − |r12 (t)| = d − r12 (t). We define the beginning of a collision at t = 0 and the end at t = tc
through the conditions, ξn (0) = ξn (tc ) = 0.
Simulation results for the normal restitution coeﬃcient for the models A and B are shown in Fig.
4.2 (upper and lower panel respectively). As it is seen from the figure, the restitution coeﬃcient e,
defined by Eq. (4.1) becomes negative for large incident angles γ. Such unusual behavior of e at
nanoscales may be understood if we notice that the orientation of the contact plane, characterized by
the unit vector n(t) = r12 (t)/r12 (t), significantly alters during the collision, Fig. 4.1. This is quantified
by the angle α between the initial and final orientations of n(t),cos α = n(0) · n(tc ) .

(4.2)

The dependence of α on the incident angle γ is shown in Fig. 4.3. If α is large, the normal restitution
coeﬃcients becomes negative, Fig. 4.2.

37

restitution coefficient

1
0.5
0
-0.5
-1
0 0.2 0.4 0.6 0.8 1.0 1.2 1.4

restitution coefficient

γ
1
0.5
0
-0.5
-1
0 0.2 0.4 0.6 0.8 1.0 1.2 1.4γ
Figure 4.2: Dependence on the incident angle γ of the normal restitution coeﬃcients e, and ẽ according to the standard definitions (4.4) and modified definition (4.5). Open circles and squares arerespectively the MD results for e and ẽ, while dashed and solid lines correspond to theoretical predictions. Upper panel refers to the model A and lower panel – to the model B. Note that the coeﬃcient ẽ
is always positive.

38

0.8
0.6

α 0.4
0.2
0
0 0.2 0.4 0.6 0.8 1.0 1.2 1.4γ
Figure 4.3: Dependence of the angular displacement α = arcos [n(0) · n(tc )] of the unit normal n,
Eq. (4.2), on the incident angle γ. Open squares and circles are the MD results for the model
R t A and B,
respectively. Solid and broken lines are the corresponding theoretical predictions, α = 0 c Ω(t)dt (seetext for detail).

4.3 Modified definition of e
To analyze this eﬀect, consider the relative velocity of particles at their contact,g = v12 +dd[n × ω12 ] = −ξ̇n n + r12 ṅ + [n × ω12 ] ,
2
2

(4.3)

where ω12 ≡ ω1 + ω2 and we use v12 = ṙ12 with r12 = n(d − ξn ). In the standard definition ofe and theoretical studies of an oblique impact [167], n is taken at the collision instant, that is, itsreorientation during the impact is ignored. In experiments, the normal n is also determined only
once, at the beginning of an impact [168, 169], n = n(0). Neglecting angular velocities (note that
ω1/2 (0) = 0) we find for the restitution coeﬃcient:e=−g(tc ) · n(0)ξ̇n (tc )d sin α α̇=cos α −,
g(0) · n(0)V cos γξ̇n (0)(4.4)

where we take into account that g(0) · n(0) = −ξ̇n (0) = −V cos γ < 0, that r12 (tc ) = d and thatξ̇n (tc ) < 0. For head-on collisions, when γ → 0 and α → 0 (see Fig. 4.3) the second term in the r.h.s.
of Eq. (4.4) is negligible and e is positive. For oblique impacts γ and α are large and the second term
prevails, yielding a negative e.
Hence, the negative restitution coeﬃcient is a consequence of a significant reorientation of a contact plane during a collision. For hard particles with a small collision duration the reorientation of nis small and may be neglected [167]; this usually holds true for macroscopic bodies. Nanoclusters,
however, are very soft particles with small Young’s modulus [161]. The duration of their impact tc is
relatively large and the reorientation of the contact plane is significant.
As it follows from the Eq. (4.4), the standard definition of e characterizes not only the normal
motion along n(t) (the first term in the r.h.s. of Eq. (4.4)), but also the change of the normal n(t) (thesecond term in the r.h.s. of Eq. (4.4)). Therefore, it is worth to define the restitution coeﬃcient, which
describes pure normal motion. The respective modification of the standard definition reads:ẽ = −

ξ̇n (tc )g(tc ) · n(tc ).
=g(0) · n(0)ξ̇n (0)39(4.5)

Here we use Eq. (4.3) for t = tc and take into account that ṅ · n = 0 for a unit vector n. Note, thatthe modified restitution coeﬃcient ẽ is always positive, Fig. 4.2. It can be also seen from Fig. 4.2 thatthe magnitude of ẽ for an oblique impact (for large γ) is significantly larger than that for a head-on
collision. In what follows we explain the observed behaviors of e and ẽ using a simple theoreticalmodel, based on continuum mechanics approach.

4.4 Theory of an oblique impact
Consider a non-inertial frame, rotating with the angular velocity Ω, perpendicular n, so that ṅ =
Ω × n. To compute the normal force acting between two nanoclusters we apply the impact theory for
macroscopic viscoelastic adhesive spheres [86,87]. It contains the JKR force [81], which accounts forelastic interactions via the Herzian force F H and for adhesive interactions via the Boussinesq forceF B,r
6πσ 3/2
4a3
−
a .
(4.6)
FH − FB =Dd
DIt also contains the dissipative force [87],r 12a2 3 6πσ 1/2 −
a  .
(4.7)
F D = ȧη 
Dd
2
D
Here, a is the contact radius of the colliding nanoclusters, related to the normal displacement ξn as
r
4a2
8πσDa
ξn =
−
,
(4.8)
d
3
and D = (3/2)(1 − ν2 )/Y is the elastic constant with the Young modulus Y and the Poisson ratio ν.
From the independent numerical simulations we estimate Y = 88.3 /σ3LJ and ν = 0.396 for modelA, and Y = 283 GPa and ν = 0.166 for model B [161]. The surface tension σ may be expressed via
Hamaker constant AH and the equilibrium distance between atoms at the interface z0 as σ ' AH /24πz20 .
We obtain σ = 0.0246 /σ2LJ and 0.00289 N/m for the models A and B, respectively. The dissipative
constant η, which accounts for the viscoelasticity of the particles’ material [86] is used here as a fitting√
parameter. In the present simulations a good agreement is obtained by choosing η = 0.65 σLJ m/
and 1.62 fs for models A and B, respectively.
In the non-inertial frame, the inertial force must be also taken into account. Its normal component
reads [170],
F I = 2µv12 · ṅ(t) − µx12 |ṅ(t)|2 ,(4.9)
where µ = Nm/2 is the reduced mass of the nanoclusters. If we again neglect the angular velocitiesof particles in the collision (since ω1/2 (0) = 0), that is, if we assume that the two clusters at a contact
move together as a solid dumbbell, we can exploit the conservation of the angular momentum in the
form,
2
µr12
Ω = µV sin γ d ,(4.10)
2 (t) and the inertial force,
where we take into account that n · Ω = 0. This yields Ω(t) = V sin γ d/r12

FI =µV 2 d2 2
sin γ .
3
r12

(4.11)

Combining Eqs. (4.6) – (4.11) we obtain the equation of motion for ξn :µd0 2
µV 2 d2
ξ
+
F
−
F
+
F
+
sin2 γ = 0 ,
n
H
B
D
dt2
(d − ξn )340(4.12)

where d0 /dt denotes the time derivative in the non-inertial frame. Solving Eq. (4.12) for ξn (t), we
obtain ẽ as it follows from Eq. (4.5). Taking into account that α̇ = Ω(tc ) we obtain from Eq. (4.4) the
relation between the standard and modified restitution coeﬃcients,e = ẽ cos α − tan γ sin α .
(4.13)
Rt
The last equation together with the relation α = 0 c Ω(t)dt may be used to compute the standard
coeﬃcient e. The theoretical predictions for the coeﬃcients e and ẽ are shown on the upper andlower panels of Fig. 4.2 respectively. The agreement between our theory, which has only one fitting
parameter, and MD simulations is rather good. We find that the restitution coeﬃcient of H-passivatedSi nanospheres is well reproduced by our macroscopic theory for the incident speed between 20 m/s
and 2405 m/s. If, however, the speed exceeds 2500 m/s, the nanospheres melt and fuse upon collisions
and the theory fails to describe the impact.
Figure 4.4 exhibits time development of ξn . Upper and lower panels refer to the Model A and B,
respectively. Open circles are the results of the MD simulation averaged over 100 samples and 10
samples for the Model A and B, respectively. Here, the incident angles are γ = 0.05 and γ = 1.25 for
the Model A and B, respectively. The solid lines are the solutions of the equation of motion Eq.(4.12).
However, it should be noted that, in the case of the Model A, Eq.(4.12) is scaled by the LJ parameters
and its form does not change even if the length L and time τ are rescaled to λL and λτ, respectively.
Therefore, in Fig.4.4,√ we rescaled the length and time units of the theoretical curve for the Model Aas λσLJ and λσLJ m/ with λ = 0.5.
We wish to stress that our theoretical model, developed for nanoclusters, may be relevant for
oblique collisions of macroscopic bodies, provided the re-orientation of the contact plane during theimpact is not negligible. This may happen for soft cohesive particles with a low Young modulus andlarge collision time. Relevance of the theory for collisions in wet granular systems is also expected[171, 172].

4.5 Conclusion
We perform a detailed study of the oblique impact of nanoclusters by means of Molecular Dynamics
and theoretically. In simulations we use two models, a simplified one, based on the Lennard-Jones
potential with a cohesive parameter and a realistic model for nanoclusters with covalently bonded
atoms. We detect unexpected behavior of the normal restitution coeﬃcient e, which becomes negative
for large incident angles and explain this eﬀect by the reorientation of the contact plane in the courseof collision. We propose a modified definition of the restitution coeﬃcient, ẽ, which describes only
the normal motion of particles, independently of their relative reorientation, and is always positive.
A simple relation between e and ẽ, that may be helpful for experiments is reported. We develop atheoretical model for an oblique impact, based on the continuum mechanics description of colliding
particles, and demonstrate that theoretical predictions agree well with simulation results. Hence,
we conclude that the macroscopic concepts of elasticity, surface tension and bulk viscosity are wellapplicable for nano-objects of a few hundreds atoms.

41

3
2
1
0

0

4

8

12

time
2
1.5
1
0.5
0

0 2 4 6 8 10 12 14

time [fs]Figure 4.4: Time development of ξn . Upper and lower panels refer to the Model A and B, respectively.
Open circles are the results of the MD simulation averaged over 100 samples and 10 samples for the
Model A and B, respectively. Here, the incident angles are γ = 0.05 and γ = 1.25 for the Model A
and B, respectively. The solid lines are the theoretical predictions, however,√ we rescaled the lengthand time units of the theoretical curve for the Model A as λσLJ and λσLJ m/ with λ = 0.5.

42Chapter 5

Summary
In Chapter 2, we find that deposited LJ nanoclusters consisting 32, 108, 255, 300, 500 and 862 atoms
exhibit a transition√ from multilayered adsorption to monolayered adsorption at the critical incident
speed, Vimp = 3.3 /m. From our simulation, we clarify that the deposited nanoclusters are perfectly
wetting on the substrate above the critical incident speed. Employing the energy conservation law, we
estimate the critical value of the incident speed at√ which the evaporation begins to occur during theimpact. The estimated critical value, Vimp = 1.7 /m, exhibits good agreement with our simulation
results. We also find that the scaled adsorption parameter is independent of the cluster size and is
proportional to the square of the incident speed. These results also exhibit good agreement with
the MD simulation results. We perform the Fourier analysis of √the boundary shape of√ the adsorbed
nanocluster and find that some modes becomes unstable for 3.0 /m < Vimp < 4.0 /m. We also
perform the MD simulation of A-C case and we find that the A cluster does not become a monolayer
film on the C surface. Thus, we conclude that not only the incident speed, but also the strengthof the interaction between the cluster and the substrate is important to form a monolayer film on
a substrate. We also investigate the temperature dependence of deposited LJ nanoclusters and findthat nanoclusters are stable if the temperature T < T C but the evaporation occurs if T > T C . The
morphology is almost independent of the temperature, and the binding energy and the orientationalbond order parameter Q6 increases as the temperature increases.
In Chapter 3, we perform the molecular dynamics simulation of the graphene sheet induced bya collision with an argon nanocluster, and the bending formation and the propagation of transverse
deflection wave are observed. We find that the linear theory of the elasticity well explains the time
evolution of the deflection of the graphene. In addition, we conclude that the actual thickness is muchthinner than the diameter of a carbon atom from the analysis of the motion of the graphene. We also
analyze the time evolution of the temperature profile, and find that the analysis based on the leastdissipation principle reproduces our simulation result in the early stage of the impact.
In Chapter 4, we perform a detailed study of the oblique impact of nanoclusters by means ofMolecular Dynamics and theoretically. In simulations, we use two models, a simplified one, basedon the LJ potential with a cohesive parameter and a realistic model for nanoclusters with covalently
bonded atoms. We detect unexpected behavior of the normal restitution coeﬃcient e, which becomes
negative for large incident angles and explain this eﬀect by the reorientation of the contact planein the course of collision. We propose a modified definition of the restitution coeﬃcient, ẽ, which
describes only the normal motion of particles, independently of their relative reorientation, and isalways positive. A simple relation between e and ẽ, that may be helpful for experiments is reported.
We develop a theoretical model for an oblique impact, based on the continuum mechanics descriptionof colliding particles, and demonstrate that theoretical predictions agree well with simulation results.
Hence, we conclude that the macroscopic concepts of elasticity, surface tension and bulk viscosity are
well applicable for nano-objects of a few hundreds atoms.

43Appendix A

Tersoﬀ potential
The Tersoﬀ potential is a bond-order potential used for the covalent bonding between two carbon(C) atoms or two silicon (Si) atoms [120–123]. Mota et. al. calculated the material parameters of the
Tersoﬀ potential to use the potential for silicon nitride and silicon hydride [124,125]. In this appendix,
we explain how to use the Tersoﬀ potential for the covalent bondings Si-Si, Si-H and H-H in the MD
simulation.

A.1 Explicit form of the Tersoﬀ potential
Interaction energy between atoms i and j is introduced as a sum of repulsive and attractive termshi
Vi j = fC (ri j ) fR (ri j ) + bi j fA (ri j ) ,(A.1)
where ri j is the distance between atoms i and j. In Eq.(A.1), fR and fA represent the repulsive andattractive interactions in the functional forms of Morse potentialfR (ri j ) = Aαβ e−λαβ ri j ,−µαβ ri jfA (ri j ) = −Bαβ e

(A.2)
,

(A.3)

and fC is a cutoﬀ function which decrease smoothly to zero between two distances Rαβ and S αβ1,ri j < Rαβ






1
π(ri j −Rαβ )(A.4)
fC (ri j ) = , Rαβ < ri j < S αβ
2 1 + cos S αβ −Rαβ




0 ,
S αβ < ri j
where the subscripts α and β represent the species of atoms i and j, respectively.
In Eq.(A.1), bi j represents the three body interaction of a triplet of three atoms i, j, and k which
are arranged as Fig.A.1. The functional form of bi j is− 1bi j = χαβ 1 + βnαα ζinjα 2nα ,X
ζi j =fC (rik )g(θi jk ) ,k,i, jcαg(θi jk ) = 1 +
dα

!2
−

c2α
,
dα2 + (hα − cos θi jk )2

(A.5)
(A.6)

(A.7)

where θi jk is the angle between i j and ik bonds, and χαβ is defined as
χSiSi = χHH = 1.0 ,(A.8)

χSiH = 0.78 .

(A.9)

44

k
θijk
j

i

Figure A.1: The configuration of three atoms i, j, and k.
Table A.1: The list of material parameters for the Tersoﬀ potential.

A (eV)
B (eV)
R (Å)
S (Å)λ (Å−1 )µ (Å−1 )
β
n
c
d
h

Si

H

1.8308 × 103

86.712043.53100.80
1.00
3.78791.98004.0000
1.00
0.00
1.00
1.00

4.7118 × 102

2.70
3.00
2.47991.73221.1000 × 10−6
7.8734 × 10−1
1.0039 × 105
1.6217 × 101−5.9825 × 10−1

The parameters in Eqs.(A.2)-(A.4) are defined asp
Aαβ =Aα Aβ ,pBαβ =Bα Bβ ,pRαβ =Rα Rβ ,p
S αβ =S αS β ,λα + λβλαβ =,2µα + µβµαβ =,2(A.10)(A.11)(A.12)(A.13)(A.14)(A.15)(A.16)

where the material parameters Aα , Bα , Rα , S α , λα , µα , βα , nα , cα , dα , and hα for silicon and hydrogenatoms (α = Si, H) are listed in Table A.1.

A.2 Forces between interacting three atoms
The bond-order coeﬃcient bi j is generally asymmetry about i and j, i.e. bi j , b ji . However, Tersoﬀ
assumed in his article that bi j and b ji can be replaced byb∗i j =

bi j + b ji.
2

(A.17)Let us introduce the positions of atoms i, j, and k as ri , r j , and rk , respectively. The forces acted45

on the atoms i, j, and k are! #∂b∗i jri j∂ fR∗ ∂ fA= −++ bi jfC− fC fA
,
∂ri j
∂ri j
∂ri j
ri j
∂ri
"
! #
∂b∗i j ∂ fCri j∂ fR∗∗ ∂ fA=fR + bi j fA+
+ bi j− fC fA
,
fC
∂ri j
∂ri j
∂ri j
ri j
∂r j
∂b∗i j= − fC fA
,
∂rk
"

∂Vi j
Fi = −∂ri
Fj = −∂Vi j∂r j

Fk = −∂Vi j∂rk

fR + b∗i j fA

 ∂ fC(A.18)(A.19)(A.20)

respectively, where ri j ≡ ri − r j . In Eqs.(A.18)-(A.20), the derivatives of fR (ri j ) and fA (ri j ) are∂ fR (ri j )∂ri j∂ fA (ri j )∂ri j= −λαβ Ae−λαβ ri j ,(A.21)= µαβ Be−µαβ ri j ,(A.22)

respectively. The derivatives of the cutoﬀ function fC (ri j ) is



0,



∂ fC (ri j ) π(r −R )=,− 2(S π−R ) sin S αβi j −Rαβ

αβ
αβ
αβ

∂ri j



0 .
From Eq.(A.17),Rαβ < ri j < S αβ

(A.23)S αβ < ri j!
1 ∂bi j ∂b ji
+
,
2 ∂ri
∂ri

(A.24)

bi j βnαα ζinjα −1 ∂ζi j
∂bi j
=−
.
∂ri
2(1 + βnαα ζinjα ) ∂ri(A.25)∂b∗i j∂ri

where

ri j < Rαβ=The derivatives of ζi j are
∂ζi j
∂ri

X"

#
∂g(θi jk )∂ fC (rik ) rik=g(θi jk )+ fC (rik ),∂rik rik∂rik,i, j

∂ζi j
∂r j

=

∂ζi j
∂rk

= −g(θi jk )X

fC (rik )k,i, j∂g(θi jk ),∂r j(A.26)(A.27)∂g(θi jk )∂ fC (rik ) rik+ fC (rik ),∂rik rik∂rk(A.28)

and we can find that∂g(θi jk )dg(θi jk ) ∂ cos θi jk
2c2α (cos θi jk − hα ) ∂ cos θi jk
=
=h
.
i2 ∂r
∂ri
d(cos θi jk ) ∂riidα2 + (hα − cos θi jk )2
Becausecos θi jk =2 − r2ri2j + rikjk

2ri j rik46,

(A.29)(A.30)

Eq.(A.29) can be calculated by using the following derivatives!
!
∂ cos θi jk
cos θi jk ri jcos θi jk rik
1
1
=
−
+
−
,
∂ri
rik
ri j
ri j
ri j
rik
rik
cos θi jk ri j
∂ cos θi jk
1 rik
=
−
,
∂r j
ri j ri j ri j rik∂ cos θi jk
cos θi jk rik
1 ri j=−.
∂rk
rik rik rik ri j

47

(A.31)(A.32)(A.33)

Appendix B

Brenner potential
The Brenner potential is also a bond-order potential used for the covalent bonding between two carbon(C) atoms [126–129]. In this appendix, we follow Ref. [126] and explain how to use the Brennerpotential for the covalent bondings C-C in the MD simulation.

B.1 Explicit form of the Brenner potential
Interaction energy between atoms i and j is introduced as a sum of repulsive and attractive termshi
Vi j = fC (ri j ) fR (ri j ) + bi j fA (ri j ) ,(B.1)
where ri j is the distance between atoms i and j. In Eq.(B.1), fR and fA represent the repulsive andattractive interactions in the functional forms of Morse potential
A −λ √2β(ri j −re )e,
β−1βA −λ √2/β(ri j −re )efA (ri j ) = −,β−1
fR (ri j ) =and fC is a cutoﬀ function which decrease smoothly to zero between two distances R and S
1,ri j < R π(ri j −R) o1 n
fC (ri j ) = , R < ri j < S
2 1 + cos
S −R



0.
S < ri j

(B.2)
(B.3)

(B.4)

In Eq.(B.1), bi j represents the three body interaction of a triplet of three atoms i, j, and k which
are also arranged as Fig.A.1. The functional form of bi j is

−δ


X
bi j = 1 +fC (rik )g(θi jk ) ,k,i, j"# c 2c2g(θi jk ) = α0 1 +− 2,
dd + (1 + cos θi jk )2

(B.5)
(B.6)

where θi jk is the angle between i j and ik bonds. The material parameters are listed in Table B.1, and
we also replace bi j and b ji by
bi j + b ji
b∗i j =
.
(B.7)
2

48

B.2 Forces between interacting three atoms
Let us introduce the positions of atoms i, j, and k as ri , r j , and rk , respectively. The forces acted on
the atoms i, j, and k are∂Vi j
Fi = −∂ri
Fj = −∂Vi j∂r j

Fk = −

∂Vi j
∂rk

"

! #
∂b∗i jri j∂ fR∗ ∂ fA= −++ bi jfC− fC fA
,
∂ri j
∂ri j
∂ri j
ri j
∂ri
"
! #
∂b∗i j ∂ fCri j∂ fR∗∗ ∂ fAfR + bi j fA=+
+ bi j− fC fA
,
fC
∂ri j
∂ri j
∂ri j
ri j
∂r j
∂b∗i j= − fC fA,∂rk
fR + b∗i j fA

 ∂ fC(B.8)(B.9)(B.10)

respectively, where ri j ≡ ri − r j . In Eqs.(B.8)-(B.10), the derivatives of fR (ri j ) and fA (ri j ) are∂ fR (ri j )∂ri j∂ fA (ri j )∂ri j= −Be−λ= Be−λ√2β(ri j −re )√2/β(ri j −re ),(B.11),(B.12)

pwhere we introduced B ≡ λ 2βA/ (β − 1). The derivatives of the cutoﬀ function fC (ri j ) is0,
ri j < R∂ fC (ri j ) 

π
r
−R
(
)
i
j
=
− 2(Sπ−R) sin S −R , R < ri j < S


∂ri j


0 .
S < ri j(B.13)!
1 ∂bi j ∂b ji
=
+
,
∂ri
2 ∂ri
∂ri

(B.14)∂bi j
bi j δ ∂ζi j=−.
∂ri
1 + ζi j ∂ri(B.15)

From Eq.(B.7),∂b∗i jwhere

In Eq.(B.15), we introduced ζi j ≡Pk,i, j fC (rik )g(θi jk ). The derivatives of ζi j are

∂ζi j
∂ri

X"

#
∂g(θi jk )∂ fC (rik ) rik=g(θi jk )+ fC (rik ),∂rik rik∂rik,i, j

∂ζi j
∂r j

=

∂ζi j
∂rk

= −g(θi jk )X

fC (rik )k,i, j∂g(θi jk ),∂r j(B.16)(B.17)∂g(θi jk )∂ fC (rik ) rik+ fC (rik ),∂rik rik∂rk(B.18)

and we can find thatdg(θi jk ) ∂ cos θi jk
2α0 c2 (1 + cos θi jk ) ∂ cos θi jk∂g(θi jk )
=
=h
,
i2 ∂r
∂ri
d(cos θi jk ) ∂riid2 + (1 + cos θ )2
i jk

which can be calculated by using Eqs.(A.31)-(A.33).

49(B.19)

Table B.1: The list of material parameters for the Brenner potential.
A(eV)re (Å)R(Å)
S (Å)
λ(Å−1 )

β
δ
α0
c
d

6.325
1.315
1.7
2.0
1.5

50

1.29
0.80469
0.011304
19.0
2.5

Appendix C

Spin eﬀect on the restitution coeﬃcient of
nanocluster
In this appendix, we derive the equation of motion Eq.(4.12) under spin-less approximation. We also
explain how to calculate the angular displacement α and derive Eqs.(4.4), (4.5) and (4.13). Finally,
we discuss the spin eﬀect on the restitution coeﬃcient of nanocluster.

C.1 Spin-less approximation
In the inertial frame, the angular momentum of the system is divided into two parts as
L = µr12 (n × v12 ) + I1 ω1 + I2 ω2≡ Lorbital + L spin ,(C.1)

where we introduce the orbital part Lorbital ≡ µr12 (n × v12 ) and the spin part L spin ≡ I1 ω1 + I2 ω2 ofthe angular momentum. In Eq.(C.1), µ and n are the reduced mass of cluster 1 and 2, and the normal
vector projected from the center of mass (CM) of cluster 1 to the CM of cluster 2, respectively. r12
and v12 are the distance between the CM of cluster 1 and the CM of cluster 2, and the translational
relative velocity, respectively. Iα and ωα are the moment of inertia and the angular velocity of cluster
α. We assume that the rotation of cluster 1 and 2 occurs in the xy plane. This means that the angularmomentum of the system only has the z-component. Thus, L = (0, 0, Lz ), Lorbital = (0, 0, Lzorbital ) andspin
L spin = (0, 0, Lz ). If we introduce the angular velocity of n as Ω = (0, 0, θ̇), we can find that[n × v12 ]z = r12 θ̇ .

(C.2)

2
Lzorbital = µr12
θ̇ .

(C.3)

Therefore, the orbital part of L becomes

Before the collision, the cluster 1 moves to the cluster 2 with the incident velocity v12 (0) = (0, −V, 0),
the angular velocities ω1 = ω2 = 0, and the incident angle γ. Thus, the angular momentum of thesystem before the collision is
L = µr12 [n(0) × v12 (0)]= −µVd sin γez ,(C.4)

here ez is the unit vector projected to the z-axis. From the conservation law of the angular momentumof the system,
Lz = −µVd sin γ .
(C.5)
51

Therefore, from Eqs.(C.1), (C.3) and (C.5), the angular velocity θ̇ becomesspin

Lz − Lz
2
µr12

θ̇ =

spin

= −

µVd sin γ + Lz
2
µr12

.

(C.6)

In the non-inertial frame, the equation of motion isµd0 0
v = Fn + F I ,
dt n

(C.7)

where d0 /dt represents the time derivative in the non-inertial frame, and v0n = d0 r12 /dt. The normal
force Fn acts between the CM of cluster 1 and the CM of cluster 2, and the inertial force for the normal
motion F I is given by
F I = 2µv12 · ṅ(t) − µr12 |ṅ(t)|2 .
(C.8)
Because |ṅ| = |Ω × n| = θ̇ andv12 · ṅ = v12 · [Ω × n]
= Ω · [n × v12 ]= r12 θ̇2 ,(C.9)

F I = 2µr12 θ̇2 − µr12 θ̇2 = µr12 θ̇2 .

(C.10)

d0 0
v = Fn + µr12 θ̇2 .
dt n(C.11)

we findThus, Eq. (C.7) becomesµspin

If we neglect Lz
given byin Eq. (C.6), which we call spin-less approximation, the angular velocity of n isVdθ̇ = − 2 sin γ .
r12(C.12)

Because the normal compression is defined as ξn = d − r12 , i.e. d0 v0n /dt = −d02 ξn /dt2 , and we adopt
the model for dissipative particles, i.e. Fn = F H − F B + F D , Eq.(C.11) gives the equation of motionof ξnµV 2 d2
d02
sin2 γ = 0 ,(C.13)µ 2 ξn + F H − F B + F D +dt(d − ξn )3which is Eq.(4.12) in Chap.4.

C.2 Angular displacement α
Because of ṅ = Ω × n, we can find the time evolution of n(t) = (n x , ny , nz ) by solving the following
equations.
Vdsin γny ,2
r12
Vd= θ̇n x = − 2 sin γn x ,
r12
= 0.

ṅ x = −θ̇ny =(C.14)ṅy(C.15)ṅz

52(C.16)

Here, the time evolution of r12 = d − ξn is found by solving Eq.(C.13). Therefore, the angulardisplacement α
cos α = n(0) · n(tc )(C.17)
can be calculated and the results are displayed in Fig.4.3.

C.3 Derivations of Eqs.(4.4), (4.5) and (4.13)
The inner products of g(0), g(tc ) and n(0) are
g(0) · n(0) = −ξ̇n (0) ,(C.18)g(tc ) · n(0) = −ξ̇n (tc ) (n(tc ) · n(0)) + d (ṅ(tc ) · n(0)) .

(C.19)

The angular displacement at time t iscos α(t) = n(0) · n(t) ,(C.20)−α̇(t) sin α(t) = n(0) · ṅ(t) .

(C.21)and the time derivative of α(t) is

Then, the restitution coeﬃcient in the inertial frame isg(tc ) · n(0)g(0) · n(0)ξ̇n (tc )d sin αα̇(tc )= −cos α −,
ξ̇n (0)ξ̇n (0)e = −(C.22)

where we omit the argument of α(tc ) ≡ α. Because ξ̇n (tc )/ξ̇n (0) is always negative and ξ̇n (0) = V cos γ,
Eq.(C.22) is rewritten as
d sin αα̇(tc )ξ̇n (tc )e=cos α −,
(C.23)
V cos γξ̇n (0)which is Eq.(4.4) in Chap.4. In the same way,ẽ = −g(tc ) · n(tc )ξ̇n (tc )=,g(0) · n(0)ξ̇n (0)(C.24)

2 ,which is Eq.(4.5) in Chap.4. Because α̇(t) = −θ̇ = Vd sin γ/r12α̇(tc ) =V sin γ.
d(C.25)

From Eqs. (C.23)-(C.25),
d sin α V sin γ
V cos γ d= ẽ cos α − tan γ sin α ,e = ẽ cos α −which is Eq.(4.13) in Chap.4.

53(C.26)

0.8
0.7
0.6
0.5
0.4
0.3
0 0.2 0.4 0.6 0.8 1.0 1.2 1.4

γ
1.0
0.9
0.8
0.7
0 0.2 0.4 0.6 0.8 1.0 1.2 1.4γ
Figure C.1: The results of the restitution coeﬃcient ẽ under the spin-less approximation (red solid
spinlines) and the results of ẽ with Lz which is measured in the MD simulation (green broken lines).
Upper and lower panels refer to the model A and B, respectively.

C.4 Spin eﬀect on the restitution coeﬃcient√
If we use the spin-less approximation, the dissipative constants are η = 0.65σ m/ and 1.62 fs formodels A and B, respectively. Figure C.1 displays the results of the restitution coeﬃcient ẽ underspin
the spin-less approximation and the results of ẽ with Lz which is measured in the MD simulation.
Upper and lower panels refer to the model A and B, respectively. When the MD simulation result
spinof Lz is introduced, ẽ for the model A is almost same with the result of spin-less approximation.
However, ẽ for the model B becomes larger than the case of spin-less approximation. In this case,the result of ẽ can be fitted with larger dissipative constant η spin = 1.68 fs. Therefore, if we use thespin-less approximation, the eﬀect of spin of nanocluster is absorbed into the dissipative constant.

54

Appendix D

Elastic modulus and surface tension of ananocluster
D.1 Elastic modulus
If an elastic body is deformed, the variation of free energy per unit volume is represented as
1
K∆F = u2ll + λ uik − δik ull
2
3

!2
,

(D.1)

where uik , K and λ are the strain tensor, the compressibility and the shear rate, respectively. In thissection, we explain how to measure Young’s modulus Y and the Poisson ratio ν of nanocluster, which
are related to K and λ by
Y =
ν =

9λK
,
3K + λ
1 3K − 2λ,2 3K + λ

(D.2)
(D.3)

respectively.

D.1.1

Isotropic compression

To measure the compressibility K of nanocluster, we consider isotropic compression. At first, we put
the center of mass of nanocluster onto the origin. If we give the displacement (u x , uy , uz ) = (ax, ay, az)
to each atom in nanocluster with the deformation rate a, the components of the strain tensor are givenbyu xx =
uyy =
uzz =

∂u x
=a,
∂x
∂uy
=a,
∂y
∂uz
=a,
∂z

(D.4)
(D.5)
(D.6)

and uik = 0 (i , k), respectively. In this case,2u2ll = u xx + uyy + uzz = 9a2 ,!21
uik − δik ull = 0 ,
3
55

(D.7)
(D.8)

where we use δll = 3. Therefore, from Eq. (D.1), the compressibility K of nanocluster is given by
K=

2∆F
.
9a2

(D.9)

If we assume the deformation of nanocluster is taken place adiabatically, the variation of freeenergy equals the variation of internal energy ∆E, which is the change of the potential energy ofnanocluster. In Chap.4, we measure ∆E in another MD simulation and find K of nanocluster fromK = 2∆E/9a2 .
Shearing deformation
To measure the shear rate λ of nanocluster, we consider shear deformation. We also put the center ofmass of nanocluster onto the origin and give the displacement (u x , uy , uz ) = (ay, 0, 0) to each atom in
nanocluster with the deformation rate a. In this case, the components of strain tensor u xy and uyx are
given by!
a
1 ∂u x ∂uy+= ,(D.10)
u xy = uyx =
2 ∂y
∂x
2
and other components are zero. Thus,
u2ll = 0 ,!21a2
uik − δik ull.
=
3
2(D.11)(D.12)

Therefore, from Eq.(D.1), the shear rate λ of nanocluster is given by
λ=

2∆F
.
a2

(D.13)In Chap.4, we also assume the deformation is adiabatic. Then, we measure ∆E in another MD
simulation and find λ of nanocluster from λ = 2∆E/a2 .
In the same way, we also give the displacements (u x , uy , uz ) = (az, 0, 0), (u x , uy , uz ) = (0, az, 0),(u x , uy , uz ) = (0, ax, 0), (u x , uy , uz ) = (0, 0, ax) and (u x , uy , uz ) = (0, 0, ay) to nanocluster, and measure
λ from λ = 2∆E/a2 . Finally, we determine λ as the averaged value of these six measurements.

D.2 Surface tension
D.2.1

Model of interacting two spheres

The surface tension σ is defined asE(leq ),
(D.14)
2
where E and leq are the interaction energy and the equilibrium distance between two nanoclusters,
respectively. The interaction energy is represented by
Z
Z
2
E=ρ
dr1
dr2 φ12 ,(D.15)σ=

V1

V2

where Vα (α = 1, 2) is the volume of nanocluster α, and φ12 is the interatomic two body potential. In
Eq.(D.15), we assume the number density of atoms in nanocluster ρ is constant.

56

20

total
repulsive
attractive

15

l
E

10

5

0

-5
1

2

3

4

5

l
Figure D.1: The red solid, the green broken, and the blue broken lines are respectively the interactionenergy E, the repulsive energy Erep , and the attractive energy Eatt as the functions of the distance l.
Here, the radius of the interacting two nanocluster is 5.23σLJ .
In the case of Model A in Chap.4, the interatomic potential φ12 is the modified LJ potential withthe cohesive parameter C = 0.2. The repulsive and attractive parts of φ12 (r) = φrep (r) + φatt (r) are
φrep (r) = 4 σ 12LJ,r
 σ 6LJ
φatt (r) = −4C,rrespectively. Thus, the interaction energy E = Erep + Eatt is represented byZ
Z
1Erep = 4ρ2 σ12dr
dr2 12 ,
1
LJ
r
V1
V2
Z
Z
1
Eatt = −4ρ2 σ6LJ C
dr1
dr2 6 .
r
V1
V2

(D.16)(D.17)(D.18)(D.19)

If we assume the interacting two nanoclusters are spheres with radii R1 and R2 , Eqs.(D.18) and (D.19)
can be integrated by the method of Ref. [173], and we find
Erep =Eatt

π2 σ12LJ C
540D

Z D+R2
D−R2"#
r + 9R1
r − 9R1
dr−,(r + R1 )9 (r − R1 )9(D.20)"#
2π2 ρ2 σ6LJ C
2R1 R2
2R1 R2D2 − (R1 + R2 )2= −++ ln 2, (D.21)3
D2 − (R1 + R2 )2 D2 − (R1 − R2 )2D − (R1 − R2 )2where D ≡ R1 +R2 +l and l is the distance between two surfaces of nanocluster (see Fig.D.1). Although
Eq.(D.20) can be integrated analytically, the result is too complicated to write this appendix. In our
MD simulation, the radius of nanocluster is R1 = R2 = 5.23σLJ .
Figure D.1 shows the results of the interaction energy E, the repulsive energy Erep and the attractive energy Eatt as the functions of the distance l, respectively. From this result, the equilibrium
distance is leq = 2.1σLJ and the surface tension σ is estimated asσ=E(leq )
' 0.54287/σ2LJ .
2(D.22)

Because Eq.(D.15) is represented by the two body interatomic potential φ12 , we can not use
Eq.(D.15) in the case of H-passivated Si nanoclusters (Model B in Chap.4) where we use the Tersoﬀ potential with the three body interaction terms.

57

D.2.2Model of compressed two spheres

If we follow the JKR theory [81], the surface tension σ is defined asE(ξn ),
2πa2

σ=(D.23)

where E(ξn ), ξn and a are the interaction energy, the normal compression and the contact radius ofcolliding two nanoclusters, respectively. If D0 /R  1, where D0 = 2 · 21/6 σLJ and R is the radius of
nanocluster, Eq.(D.23) can be reduced toσ=

AH,
24πD20(D.24)

where AH is the Hamaker constant.
In the case of Model A in Chap.4, we can find
AH = 4π2Cρ2 σ6LJ ,(D.25)and calculate the surface tension σ analytically. However, in the case of Model B in Chap.4, we can
not calculate AH analytically. Therefore, we measure E(ξn ) and a by another MD simulation and find
σ of H-passivated Si nanocluster from Eq.(D.23).

58

Bibliography59[1] M. Terrones et. al. Nano Today 5 (2010), 351.
[2] G.D. Scholes and G. Rumbles, Nature Materials. 5 (2006), 683.
[3] H. Dai, Acc. Chem. Res. 35 (2002), 1035.
[4] T. L. Hill, Thermodynamics of Small Systems, Parts I and II (Benjamin, Amsterdam, 1964).
[5] B.M. Smirnov and R.S. Berry, Phase Transitions of Simple Systems (Springer-Verlag, Berlin and
Heidelberg, 2007).
[6] F. Baletto and R. Ferrando, Rev. Mod. Phys. 77 (2005), 371.
[7] R.S. Berry and B.M. Smirnov, Physics Uspekhi 48 (2005), 345.
[8] A. Rytkönen, S. Valkealahti and M. Manninen, J. Chem. Phys. 106 (1997), 1888.
[9] A. Rytkönen, S. Valkealahti and M. Manninen, J. Chem. Phys. 108 (1998), 5826.
[10] C.L. Cleveland, U. Landman, T.G. Schaaﬀ, M.N. Shafigullin, P.W. Stephens and R.L. Whetten,
Phys. Rev. Lett. 79 (1997), 1873.
[11] H.-S. Nam, N.M. Hwang, B.D. Yu and J.-K. Yoon, Phys. Rev. Lett. 30 (2002), 275502.
[12] K. Koga, T. Ikeshoji and K. Sugawara, Phys. Rev. Lett. 92 (2004), 115507.
[13] C.-L. Kuo and P. Clancy, J. Phys. Chem. B 109 (2005), 13743.
[14] M.Á G.-Pinilla et. al. J. Phys. Chem. C 112 (2008), 13492.
[15] E.G. Noya, J.P.K. Doye and F. Calvo, Phys. Rev. B 73 (2006), 125407.
[16] F. Baletto, C. Mottet and R. Ferrando, Chem. Phys. Lett. 354 (2002), 82.
[17] Y. Qi et. al. J. Chem. Phys. 115 (2001), 385.
[18] J.P.K. Doye and D.J. Wales, J. Chem. Phys. 102 (1995), 9659.
[19] J.P.K. Doye, D.J. Wales and R.S. Berry, J. Chem. Phys. 103 (1995), 4234.
[20] J.P.K. Doye and D.J. Wales, Phys. Rev. Lett. 80 (1998), 1357.
[21] J.P.K. Doye, D.J. Wales and M.A. Miller, J. Chem. Phys. 109 (1998), 8143.
[22] J.P.K. Doye and F. Calvo, Phys. Rev. Lett. 86 (2001), 3570.
[23] J.P.K. Doye, Phys. Rev. Lett. 88 (2002), 238701.
[24] J.P.K. Doye and D.J. Wales, J. Chem. Phys. 116 (2002), 3777.
[25] J.P.K. Doye and F. Calvo, J. Chem. Phys. 116 (2002), 8307.
[26] F. Calvo, J.P.K. Doye and D.J. Wales, Chem. Phys. Lett. 366 (2002), 176.
[27] F. Calvo and J.P.K. Doye, Phys. Rev. B 69 (2004), 125414.
[28] J.P.K. Doye and C.P. Massen, J. Chem. Phys. 122 (2005), 084105.
[29] E.G. Noya and J.P.K. Doye, J. Chem. Phys. 124 (2006), 104503.

60

[30] N.V. Krainyukova, Eur. Phys. J. D 43 (2007), 45.
[31] W. Lechner and C. Dellago, J. Chem. Phys. 129 (2008), 114707.
[32] D. J. Wales, Energy Landscapes with Applications to Clusters, Biomolecules and Glasses (Cambridge University, Cambridge, England, 2003).
[33] Y. Imry, Phys. Rev. B 21 (1980), 2042.
[34] P. Labastie and R.L. Whetten, Phys. Rev. Lett. 65 (1990), 1567.
[35] M. Schmidt, R. Kusche, B. von Issendorﬀ and H. Haberland, Nature 393 (1998), 238.
[36] M. Schmidt, R. Kusche, T. Hippler, J. Donges, W. Kronmüller, B. von Issendorﬀ and H. Haberland, Phys. Rev. Lett. 86 (2001), 1191.
[37] W.A. de Heer, Rev. Mod. Phys. 65 (1993), 611.
[38] C. Binns, Surf. Sci. Rep. 44 (2001), 1.
[39] K. Wegner, P. Piseri, H.V. Tafreshi and P. Milani, J. Phys. D: Appl. Phys. 39 (2006), R439.
[40] O.G. Danylchenko, S.I. Kovalenko and V.N. Samovarov, Low Temp. Phys. 30 (2004), 743.
[41] W. Polak and A. Patrykiejew, Phys. Rev. B 67 (2003), 115402.
[42] N. Nishio, J. Koga, T. Yamaguchi and F. Yonezawa, J. Phys. Soc. Jpn. 73 (2004), 627.
[43] W. Polak, Eur. Phys. J. D 40 (2006), 231.
[44] W. Polak, Phys. Rev. E 77 (2008), 031404.
[45] V. V. Hoang and T. Odagaki, Phys. Rev. B 77 (2008), 125434.
[46] T. Bachels and R. Schäfer, Chem. Phys. Lett. 324 (2000), 365.
[47] M.O. Watanabe, N. Uchida and T. Kanayama, Phys. Rev. B 61 (2000), 7219.
[48] W.W. Gerberich et. al. J. Mech. Phys. Solids 51 (2003), 979.
[49] S.L. Girshick and J. Hafiz, J. Phys. D: Appl. Phys. 40 (2007), 2354.
[50] P. Jensen, Rev. Mod. Phys. 71 (1999), 1695.
[51] I. Yamada, H. Inokawa and T. Takagi, J. Appl. Phys. 56 (1984), 2746.
[52] I. Yamada, G. H. Takaoka, H. Usui, F. Satoh, Y. Itoh, K. Yamashita, S. Kitamoto, Y. Namba, Y.
Hashimoto, Y. Maeyama and K. Machida, Nucl. Instrum. Methods Phys. Res. B 59 (1991), 216.
[53] T. Takagi, Z. Phys. D 3 (1986), 271.
[54] X. Kang, L. T. Wille, H. Dreyssé and J. Eugène, Nucl. Instrum. Methods Phys. Res. B 122(1997), 339.
[55] K. Nordlund, T. T. Järvi, K. Meinander and J. Samela, Appl. Phys. A 91 (2008), 561.
[56] C.L. Cleveland and U. Landman, Science 257 (1992), 355.
[57] Y. Guoqing et. al. Nucl. Inst. Meth. Phys. Res. B 135 (1998), 382.

61

[58] K. Koga and K. Sugawara, Surf. Sci. 529 (2003), 23.
[59] K-H. Müller, J. Appl. Phys. 61 (1987), 2516.
[60] Y. Yamamura, Nucl. Instrum. Methods Phys. Res. B 45 (1990), 707.
[61] H. Haberland, Z. Insepov and M. Moseler, Phys. Rev. B 51 (1995), 11061.
[62] Y. Yamaguchi and J. Gspann, Phys. Rev. B 66 (2002), 155408.
[63] R. Krämer, Y. Yamaguchi and J. Gspann, Surf. Interface. Anal. 36 (2004), 148.
[64] A. Terasaki, J. Phys. Chem. A 111 (2007), 32.
[65] S. Bouneau, A. Brunelle, S. Della-Negra, J. Depauw, D. Jacquet, Y. LeBeyec, M. Pautrat, M.
Fallavier, J. C. Poizat and H. H. Andersen, Phys. Rev. B 65 (2002), 144106.
[66] J. Samela, J. Kotakoshi, K. Nordlund and J. Keinonen, Nucl. Instrum. Methods Phys. Res. B 239(2005), 331.
[67] M. Moseler, O. Rattunde, J. Nordiek and H. Haberland, Nucl. Instrum. Methods Phys. Res. B
164-165 (2000), 522.
[68] A. Tomsic, P. U. Andersson, N. Markovic, W. Piskorz, M. Svanberg and J. B. C. Pettersson, J.
Chem. Phys. 115 (2001), 10509.
[69] T. T. Järvi, A. Kuronen, K. Meinander, K. Nordlund and K. Albe, Phys. Rev. B 75 (2007),
115422.
[70] A. Awasthi, S. C. Hendy, P. Zoontjens and S. A. Brown, Phys. Rev. Lett. 97 (2006), 186103.
[71] A. Awasthi, S. C. Hendy, P. Zoontjens, S. A. Brown and F. Natali, Phys. Rev. B 76 (2007),
115437.
[72] M. Suri and T. Dumitricǎ, Phys. Rev. B 78 (2008), 081405(R).
[73] A. Castellanos, Adv. Phys. 54 (2005), 263.
[74] J. Tomas, Chem. Eng. Sci. 62 (2007), 1997.
[75] J. Tomas, Chem. Eng. Sci. 62 (2007), 5925.
[76] R. Tykhoniuk, J. Tomas, S. Luding, M. Kappl, L. Heim and H.-J. Butt, Chem. Eng. Sci. 62(2007), 2843.
[77] A. Bertram and J. Tomas, Micro-Macro-Interactions in Structured Media and Particle Systems(Springer-Verlag, Berlin and Heidelberg, 2008).
[78] H. Hertz, J. reine und angewandte Mathematik, 92 (1882), 156.
[79] K.L. Johnson, Contact mechanics (Cambridge University Press, The Edinburgh Building, Cambridge CB2 2RU, UK, 1985).
[80] D.A. Hills, D. Nowell and A. Sackfield, Mechanics of Elastic Contact (Butterworth-HeinemannLtd, Linacre House, Jordan Hill, Oxford OX2 8DP, UK, 1993).
[81] K. L. Johnson, K. Kendall, and A. D. Roberts, Proc. R. Soc. London 324 (1971), 301.

62

[82] N. V. Brilliantov and T. Pöschel, Kinetic Theory of Granular Gases (Oxford University Press,
Oxford, 2004).
[83] A. Domenech, Int. J. Mech. Sci. 50 (2008), 752.
[84] F.G. Bridges and A. Hatzes and D.N.C. Lin, Nature (London) 309 (1984), 333.
[85] G. Kuwabara and K. Kono, Jpn. J. Appl. Phys. 26 (1987), 1230
[86] N. V. Brilliantov, F. Spahn, J. M. Hertzsch, and T. Pöschel, Phys. Rev. E 53 (1996), 5382.
[87] N. V. Brilliantov, N. Albers, F. Spahn, and T. Pöschel, Phys. Rev. E 76 (2007), 051302.
[88] M.Y. Louge and M.E. Adams, Phys. Rev. E 65 (2002), 021303.
[89] H. Kuninaka and H. Hayakawa, Phys. Rev. Lett. 93 (2004), 154301.
[90] H. Hayakawa and H. Kuninaka, Phase Transit. 77 (2004), 889.
[91] M. Kalweit and D. Drikakis, J. Comput. Theor. Nanosci. 1 (2004), 367.
[92] M. Kalweit and D. Drikakis, Phys. Rev. B 74 (2006), 235415.
[93] H. Kuninaka and H. Hayakawa, Phys. Rev. E 79 (2009), 031309.
[94] A. Awasthi, S.C. Hendy, P. Zoontjens and S.A. Brown, Phys. Rev. Lett. 97 (2006), 186103.
[95] A.I. Ayesh, S.A. Brown, A. Awasthi, S.C. Hendy, P.Y. Convers, and K. Nichol, Phys. Rev. B 81
(2010), 195422.
[96] C. Bostedt et. al. Phys. Rev. Lett. 100 (2008), 133401.
[97] A. K. Geim and K. S. Novoselov, Nature materials 6 (2007), 183.
[98] M. H. Gass, U. Bangert, A. L. Bleloch, P. Wang, R. R. Nair, and A. K. Geim, Nature nanotechnology 3 (2008), 676.
[99] D. Garcia-Sanchez, A. M. van der Zande, A. S. Paulo, B. Lassagne, P. L. McEuen, and A. Bachtold, Nano letters 8 (2008), 1399.
[100] J. C. Meyer, A. K. Geim, M. I. Katsnelson, K. S. Novoselov, T. J. Booth, and S. Roth, Nature
446 (2007), 60.
[101] J. S. Bunch, A. M. van der Zande, S. S. Verbridge, lan W. Frank, D. M. Tanenbaum, J. M.
Parpia, H. G. Craighead, and P. L. McEuen, Science 315 (2007), 490.
[102] Y. Huang, J. Wu, and K. C. Hwang, Phys. Rev. B 74 (2006), 245413.
[103] C. D. Reddy, S. Rajendran, and K. M. Liew, Nanotechnology 17 (2006), 864.
[104] Q. Lu, M. Arroyo, and R. Huang, J. Phys. D : Appl. Phys. 42 (2009), 102002.
[105] J. W. Jiang, J. S. Wang, and B. Li, Phys. Rev. B 80 (2009), 113405.
[106] P. Liu and Y.W. Zhang, Appl. Phys. Lett. 94 (2009), 231912.
[107] J. Zhou and R. Huang, J. Mech. Phys. Solids 56 (2008), 1609.
[108] B. I. Yakobson, C. J. Brabec, and J. Bernholc, Phys. Rev. Lett. 76 (1996), 2511.
63

[109] C.Q. Ru, Phys. Rev. B 62 (2000), 9973.
[110] M. Arroyo and T. Belytschko, Phys. Rev. B 69 (2004), 115415.
[111] C. Tang, W. Guo and C. Chen, Phys. Rev. Lett. 100 (2008), 175501.
[112] V. B.-Koutecký, P. Fantucci and J. Koutecký, Chem. Rev. 91 (1991), 1035.
[113] I. Boustani, W. Pewestorf, P. Fantucci, V. B.-Koutecký and J. Koutecký, Phys. Rev. B 35 (1987),
9437.
[114] P. Hohenberg and W. Kohn, Phys. Rev. 136 (1964), B864.
[115] W. Kohn and L.J. Sham, Phys. Rev. 140 (1965), A1133.
[116] E. Aprà, F. Baletto, R. Ferrando and A. Fortunelli, Phys. Rev. Lett. 93 (2004), 065502.
[117] P. Nava, M. Sierka and R. Ahlrichs, Phys. Chem. Chem. Phys. 5 (2003), 3372.
[118] D.M. Ceperley, Rev. Mod. Phys. 67 (1995), 279.
[119] L. Mitás, J.C. Grossman, I. Stich and J. Tobik, Phys. Rev. Lett. 84 (2000), 1479.
[120] J. Tersoﬀ, Phys. Rev. B 37 (1988), 6991.
[121] J. Tersoﬀ, Phys. Rev. B 38 (1988), 9902.
[122] J. Tersoﬀ, Phys. Rev. Lett. 61 (1988), 2879.
[123] J. Tersoﬀ, Phys. Rev. B 39 (1989), 5566(R).
[124] F.B. Mota, J.F. Justo and A. Fazzio, Phys. Rev. B 58 (1998), 8323.
[125] F.B. Mota, J.F. Justo and A. Fazzio, J. Appl. Phys. 86 (1999), 1843.
[126] D.W. Brenner, Phys. Rev. B 42 (1990), 9458.
[127] S.J. Stuart, A.B. Tutein and J.A. Harrison, J. Chem. Phys. 122 (2000), 6472.
[128] D.W. Brenner, O.A. Shenderova, J.A. Harrison, S.J. Stuart, B. Ni and S.B. Sinnott, J. Phys.:
Condens. Matter 14 (2002), 783.
[129] B. Ni, K.-H. Lee and S.B. Sinnott, J. Phys.: Condens. Matter 16 (2004), 7261.
[130] K. Saitoh and H. Hayakawa, Prog. Theor. Phys 122 (2009), 1081.
[131] K. Saitoh and H. Hayakawa, Phys. Rev. B 81 (2010), 115447.
[132] K. Saitoh, A. Bodrova, H. Hayakawa, and N. V. Brilliantov, Phys. Rev. Lett. 105 (2010),
238001.
[133] T. Ikeshoji, G. Torchet, M.-F. de Feraudy and K. Koga, Phys. Rev. E 63 (2001), 031101.
[134] A. A. Radzig and B. M. Smirnov, Parameters of Atoms and Atomic Ions (Énergoatomizdat,
Moscow, 1986).
[135] D. D. Do and H. D. Do, Adsorpt. Sci. Technol. 21 (2003), 389.
[136] A. Wongkoblap, D. D. Do and D. Nicholson, Phys. Chem. Chem. Phys. 10 (2008), 1106.

64

[137] P. Deltour, J-L. Barrat and P. Jensen, Phys. Rev. Lett. 78 (1997), 4597.
[138] L. J. Lewis, P. Jensen, N. Combe and J-L. Barrat, Phys. Rev. B 61 (2001), 16084.
[139] K. Yasuoka and M. Matsumoto, J. Chem. Phys. 109 (1998), 8451.
[140] J. J. Nicolas, K. E. Gubbins, W. B. Streett and D. J. Tildesly, Mol. Phys. 37 (1979), 1429.
[141] W. Christen, U. Even, T. Raz and R. D. Levine, J. Chem. Phys. 108 (1998), 10262.
[142] Q. Hou, M. Hou, L. Bardotti, B. Prével, P. Mélinon, and A. Perez, Phys. Rev. B 62 (2000),
2825.
[143] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Clarendon press, Oxford,
1987).
[144] S. V. Khare, N. C. Bartelt and T. L. Einstein, Phys. Rev. Lett. 75 (1995), 2148.
[145] P. Nozières, Solids Far From Equilibrium, Chapter 1 (Cambridge University, Cambridge, England, 1992).
[146] A. Pimpinelli and J. Villain, Physics of Crystal Growth (Cambridge University Press, Cambridge, United Kingdom, 1998).
[147] P.J. Steinhardt, D.R. Nelson and M. Ronchetti, Phys. Rev. B 28 (1983), 784.
[148] D. Cheng and D. Cao, Chem. Phys. Lett. 461 (2008), 71.
[149] M.Á. Gracia-Pinilla, E.P. Tijerina, J.A. Garcı́a, C. Fernández-Navarro, A. Tlahuice-Flores, S.
Mejı́a-Rosales, J.M. Montejano-Carrizales and M. José-Yacamán, J. Phys. Chem. C 112 (2008),
13492.
[150] E. V. Pugina, G. V. Kornich and G. Betz, Phys. Solid State 49 (2007), 580.
[151] N. Inui, K. Mochiji, and K. Moritani, Nanotechnology 19 (2008), 505501.
[152] A. A. Radzig and B. M. Smirnov, Parameters of Atoms and Atomic Ions (Énergoatomizdat,
Moscow, 1986).
[153] D. D. Do and H. D. Do, Adsorpt. Sci. Technol. 21 (2003), 389.
[154] A. Wongkoblap, D. D. Do, and D. Nicholson, Phys. Chem. Chem. Phys. 10 (2008), 1106.
[155] V. B. Shenoy, C. D. Reddy, A. Ramasubramaniam, and Y. W. Zhang, Phys. Rev. Lett. 101(2008), 245501.
[156] C. D. Reddy, A. Ramasubramaniam, V. B. Shenoy, and Y. W. Zhang, Appl. Phys. Lett. 94(2009), 101904.
[157] K. Yasuoka and M. Matsumoto, J. Chem. Phys. 109 (1998), 8451.
[158] L. D. Landau and E. M. Lifshitz, Theory of Elasticity 3rd ed (Oxford, New York, 1986).
[159] S. Timoshenko and S. W. Krieger, Theory of Plates and Shells 2nd ed. (McGraw-Hill, New
York, 1959).
[160] H. Kuninaka and H. Hayakawa, J. Phys. Soc. Japan 75 (2006), 074001.
[161] H. Kuninaka and H. Hayakawa, Prog. Theor. Phys. Suppl. 178 (2009), 157.
65

[162] H. S. Carslaw and J. C. Jaeger, Conduction of Heat in Solids 2nd ed (Oxford University Press,
Amen House, London E. C. 4, 1959).
[163] V. L. Berdichevsky, Variational Principles of Continuum Mechanics vol 1 (Springer Heidelberg
Dordrecht London New York, 2009).
[164] K. F. Riley, M. P. Hobson, and S. J. Bence, Mathematical Methods for Physics and Engineering
3rd ed (Cambridge University Press, The Edinburgh Building, Cambridge CB2 2RU, UK, 2006).
[165] A. W. Leissa, J. Sound and Vibration 31 (1973), 257.
[166] P. S. Frederiksen, J. Sound and Vibration 186 (1995), 743.
[167] T. Schwager, V. Becker, and T. Pöschel, Eur. Phys. J. E 27 (2008), 107.
[168] L. Labous, A.D. Rosato, and R.N. Dave, Phys. Rev. E 56 (1997), 5717.
[169] C.M. Sorace, M.Y. Louge, M.D. Crozier and V.H.C. Law, Mech. Res. Commun. 36 (2009),
364.
[170] L. D. Landau and E. M. Lifshitz, Mechanics (Oxford University Press, Oxford, 1965).
[171] A.A. Kantak, C.M. Hrenya and R.H. Davis, Phys. Fluids 21 (2009), 023301.
[172] C.M. Donahue, C. M. Hrenya and R. H. Davis, Phys. Rev. Lett. 105 (2010), 034501.
[173] H.C. Hamaker, Physica 10 (1937), 1058.

66