J. Plasma
J.
Plasma Fusion
Fusion Res.
Res. SERIES,SERIES,Vol.
Vol.6 6(2004)(2004)80–83
000–000

On Interatomic Potential Functions for Molecular Dynamic(MD) Simulations of Plasma-Wall Interactions
HAMAGUCHI Satoshi, OHTA Hiroaki and YAMADA Hideaki
Department of Fundamental Energy Science, Kyoto University, Uji 611-0011, Japan(Received: 9 December 2003 / Accepted: 1 March 2004)

Abstract
Control of plasma-wall interactions is a key issue in plasma material processing such as reactive ion etching(RIE) and plasma enhanced chemical vapor deposition (PECVD). For magnetic fusion devices, bulk plasma conditions
are also often influenced by interactions of the plasma with the first wall and divertor materials. Molecular dynamics(MD) simulations of plasma-wall interactions can give good insight into surface dynamics for plasma-wall interactions.
In the light of applications to classical MD simulations, we discuss interatomic potential functions that are designed
to describe covalent bonds of surface materials. Reviewing recently developed interatomic potential functions, we
examine how such potential models can represent realistic chemical bonds and how the accuracy of employedinteratomic potential functions can affect the simulation results.
Keywords:
Molecular Dynamics (MD) simulation, plasma-wall interaction, plasma-surface interaction, plasma processing,reactive ion etching

1. Introduction

simulations. Therefore in MD simulations the employedpotential model needs to be sufficiently representative of thetrue atomic interactions for the system. In this article, we
review recently developed interatomic potential functions andexamine how such potential models can represent realistic
chemical bonds and how the accuracy of employedinteratomic potential functions can affect the simulation
results.

Plasma processes such as reactive ion etching (RIE) andplasma enhanced chemical vapor deposition (PECVD)
constitute a large and essential part of semiconductor manufacturing processes today [1,2]. For such applications as wellas other industrial applications of plasmas, controlling
plasma-wall (i.e, plasma-surface) interactions is of significant
importance since the quality of processed materials is largely
determined by such interactions. For magnetic fusion devices,too, study of plasma-wall interactions has attracted muchattention of the research community as the bulk plasma
conditions are often influenced by interactions of the plasma
with the first wall and divertor materials.
Recently various investigators have examined surface
reaction mechanisms for materials exposed to beams orplasmas by using molecular dynamics (MD) simulations [39]. It has been demonstrated that such MD simulations canprovide useful, and sometimes quantitatively reliable, information on plasma-surface interaction.
For classical MD simulations, reliability of their results
is in general largely affected by the choice of interatomicpotential functions for the particles constituting the system.
To simulate plasma-surface interaction, one usually needs tohandle a large number of particles (i.e., atoms) and repeat
similar simulation runs for sufficiently many times to reduce
statistical noise, the simulation cost of which precludes theuse of more detailed ab initio (i.e., quantum mechanical)

2. Potential functions for covalent bonds
When several atoms are interacting, the force exerted oneach atom may be regarded as a multi-body force in the timescale of atomic motion since electron orbitals are often sharedby multiple atoms. If there are N atoms and the position of
the ith atom is given by ri, then the total potential energy ofthe system V may be written in cluster expansion as(1)(3)
Here vij(2) is called the two-body potential, vijk
the three-bodypotential, and so on. The functions also depend on atomic
species denoted by the subscripts. This expansion can be done
uniquely, regardless of the nature of atomic interactions. For

Corresponding author’s e-mail: hamaguch@ppl.eng.osaka-u.ac.jp180©2004 by The Japan Society of Plasma
Science and Nuclear Fusion ResearchHamaguchi S. et al., On Interatomic Potential Functions for Molecular Dynamic (MD) Simulations of Plasma-Wall Interactions

example, interactions of rare gas atoms may be well approximated by the first term only, possibly by the well-known
Lennard-Jones potential. In general, for given atoms, it isnearly impossible to construct classical interatomic potentialfunctions with high accuracy in the form of eq. (1) since it ishardly practical to take into account all possible multi-bodyinteractions. Therefore, depending on the purpose of simulations, we selectively incorporate only the important effects ofatomic interactions and represent them by relatively simple
functional forms, so that we can achieve both computational
efficiency of classical MD simulations and reasonableaccuracy of atomic dynamics.
Stillinger and Weber [3] suggested that the total potential
energy of the system may be written asF=Â v (i, j ) + Â v (i , j , k )2i< jsystem consisting of Si, F, and O atoms) and Si-Cl-O systems
are found in ref. [5].
We note that the concept of “covalent bonds” are notdirectly incorporated in the SW potential. However, the threebody potential hijk works to help atoms form effective bonds.
To see this, it is instructive to consider the exchange reaction

H + H2 Æ H2 + H ,where a single H atom is assumed to approach a hydrogenmolecule H2 in the direction of its bond axis. The parameters
for the SW potential for hydrogen are given in Table 1, which
we obtained by fitting the functions to ab initio data. [6]Figures 1 and 2 show the potential energy of the system. As
shown in Fig. 2, r12 and r23 denote the distances between theH atoms. In Fig. 1 filled points represents potential energies
obtained from ab initio quantum mechanical calculations andthe solid lines denotes the corresponding values of the SW
potential model.
It is shown in Fig. 1 that, as a single H atom approaches
the H2 molecule, i.e., r23 decreases, the potential well seen inFig. 1(a) disappears. This indicates that no bond is formed
between H and H2 and the atoms of H2 are tightly bound toeach other until the distance between H and H2 reaches about
1 Å or less. As can be seen in Figs. 1(c), if the three H atoms
are sufficiently close, no attractive force is present. The
potential energy profile as a function of r12 and r23, given by
the obtained SW potential, is shown in Fig. 2. This figure
demonstrates the formation of bonds is incorporated in theSW potential and it is controlled by the three-body potential.
There is another class of potential functions that directly
take into account the formation of covalent bonds. Based on
the general expression for bond energies derived by Abell[10], Tersoff proposed a functional form of covalent bonds[11]. The basic idea of this potential (which we call AbellTersoff potential or AT potential in short) is first to identify
bonds among atoms and then to have the strength of eachbond depend on the positions of surrounding atoms. Under
these assumptions, the total potential energy of the system
may be given by(2)3i< j <kwith v2(i, j) = vij(2)(ri, rj) being the pair interaction:

È Aij
Ê Cij ˆ Bij
Ê Dij ˆ v2 (i, j ) = Í p expÁ- q expÁ˜˜ijij
ÍÎ rijË rij - aij ¯ rij
Ë rij - aij ¯ (3)if rij < aij and v2(i, j) = 0 otherwise. Here rij = | rj – ri | is thedistance between the ith and jth atoms located at ri and rj.
The parameters, Aij, Bij, Cij, Dij, pij, qij and aij depend only on
the element types (such as silicon or oxygen) of the ith andjth atoms.
The three-body term
( 3)
v3 (i, j , k ) = vijk
( ri , rj , rk )in eq. (2) may be written asv3 (i, j , k ) = h jik (rij , rik ,q jik ) + hijk (rji , rjk ,q ijk )(4)+ hikj (rki , rkj ,q ikj )with qjik being the angle spanned by vectors rij  rj – ri andrik  rk – ri at vertex ri. The commonly used functional form
[3] for hjik is given byÊ

h jik (r, s,q ) = ÁË l jik | cos q - cos q 0jik |2a jikjÊ g jik
g kjik ˆ¥ exp Á+jk ˜Ë r - a jik s - a jik ¯(6)V=ˆ+ m jik ˜¯Â [V (r ) - B V (r )]Rijij

Aij ij(7)

iji< j(5)Table 1 Parameters of the SW potential functions, i.e., eqs. (3),(4) and (5), for hydrogen atoms. Units of energy andlength are eV and Å. Units of parameters are selectedaccordingly.

jkj
if r < ajik
and s < ajik. Otherwise hjik = 0. Here, ljik, mjik, g jik,kj
k0g jik, a jik, a jik, qjik, and a jik are parameters that depends onlyon the element types of the ith, jth, and kth atoms. Here theoriginal Stillinger-Weber (SW) functional form is modified
slightly by introducing a new parameter a jik in order toimprove parameter fitting. As can be seen from the functional
form with lijk being non negative (lijk  0), the three-bodypotential above is designed to penalize the system if atoms
are not configured to form the most stable molecules. For
example, the potential parameters for the Si-F-O (i.e., a2 bodyparameters

A HH
B HH
C HH
D HHp HHq HHa HH

281

3 bodyparameters
7.709
16.7470.73980.76631.914
1.035
2.0l HHH
mHHH
g HHH
a HHH
a HHH

0353.15
3.433
1.0
2.488

Hamaguchi S. et al., On Interatomic Potential Functions for Molecular Dynamic (MD) Simulations of Plasma-Wall Interactions

Fig. 2The potential energy of three linearly aligned H atoms
given by the SW potential. Bond formation of theclosest two H atoms is shown to be reproduced by theSW type potential due to the three-body interactionterm.

One can also combine the SW potential with the AT
potential by making the SW potential functions depend
explicitly on the conditions of surrounding atoms. For
example, one may write the total potential energy as

F=Â v (i, j; g )+Â h2

i< jFig. 1

i , j , k ( rij , rjk , q ijk ; gij , g jk ) ,ij(8)i, j ,kwhere gij represents the bond order of the possible bond
between the ith and jth atoms. (gij = 0 if there is no covalentbond between the ith and jth atoms.) This factor is typicallya function of the coordination numbers or more generally afunction of the positions of surrounding atoms. The idea is touse the two-body and three-body potentials only for the pair(ri, rj) that forms a covalent bond.
We note that the bond order gij can be determined by
how the valence electrons of the atom at each end of the bondare shared with neighboring atoms. For examples a carbon
atom has four valence electrons that may be shared withneighboring atoms. In the neighborhood of this carbon atom,
if there are more than four atoms (e.g., 5 H atoms) that canshare the electrons with this C, the same (excessive) numberof bonds (e.g., 5 C-H bonds) may be formed but the strengthof each bond is reduced accordingly. On the other hand, forexample, if there are fewer than four carbon atoms in theneighborhood of a carbon atom, some of the valence electrons
may be used to form p bonds, resulting in double or triple
bonds between carbon atoms.
Once the dependence on gij is introduced in the two-bodyterm, the separation of the three-body term hijk of eq. (8) fromthe total potential energy becomes somewhat arbitrary since

The potential energy of three linearly aligned H atoms
H-H-H. The distances between the first and second, andsecond and third H atoms are denoted by r12 and r23.
Symbol ‡ denotes data points obtained from ab initio
quantum mechanical calculations and the solid lines
are the values of eqs. (3), (4) and (5) with parameters
given in Table 1. The zero energy level is set to be thebottom of the potential well, i.e., the potential energy
for the most stable state of H2 molecule.

where V ijR denotes the repulsive potential and V ijA denotes theattractive potential. Effects of multi-body interactions areincluded in coefficient Bij. Coefficient Bij typically depends
on the numbers of atoms surrounding the atom at each endof the bond (i.e., coordination number) and the angles formedby the bond and neighboring atoms. Morse potentials [µexp(–br)] are usually used for functions V ijR and V ijA. Examplesof the AT potential for the C/H system are found in ref. [12].
382Hamaguchi S. et al., On Interatomic Potential Functions for Molecular Dynamic (MD) Simulations of Plasma-Wall Interactions

carbon atoms in the substrate material becomes significantly
large. Some details of these simulations may be found in ref.
[9] and the rest will be published elsewhere. We have also
performed MD simulations of an organic polymer, details of
which may be found in refs. [7,8].
To summarize, we have discussed interatomic potentialfunctions that are used in MD simulations for plasma (orbeam)-surface interactions. The AT potential is especiallydesign to incorporate covalent bond structures but, asdemonstrated in Fig. 2, the SW potentials can also modelbond formation with good accuracy if the three-body potential
is appropriately chosen. We can also extend the SW potentials
to explicitly incorporate bond structures as given in eq. (8).
For MD simulations of plasma-surface interactions, wherekinetic energies of impinging species are much higher thaninteratomic bond energies, physical quantities such as etchrates, deposition rates, etch product species, etc., may notdepend sensitively on details of interatomic potential functions. From various MD simulations that have been performed
to date, it seems these macroscopic quantities are mostsensitively depend on bond energies and bond lengths givenby the interatomic potentials. To confirm this, however, one
would need further analyses on the dependence of such
macroscopic data on the potential functions.

some of the three-body (and higher-order) effects are thenalready included in gij. In the case of the AT potential, thebond-angle [i.e., qijk in eq. (5) ] dependence [12] is usually
incorporated in coefficient of Bij eq. (7) whereas, in the caseof eq. (8), the bond-angle dependence is left in the “threebody” term hijk.
So far we have discussed potential functions that represents only covalent bonds. Other interatomic interactions such
as Van der Waals (VW) interactions and Coulomb interactions
can also play important roles. For example, the importanceof VW interactions is discussed in a companion paper ref.
[9].

3. Discussion
Using potential functions similar to eq. (8) we havedeveloped covalent-bond potential functions for Si/O/F/C
systems [9] and C/N/H systems [6]. Using these model
potentials together with appropriate VW interaction potentials,
we have applied MD simulations to study plasma etching anddeposition processes. An example of such simulation results
is given in Figure 3, which shows sputtering yields for SiO2
etching by F or fluorocarbon (i.e., CFx with x = 1, 2, 3)
beams. In general, sputtering yields depend on injection dose
and the results shown in Fig. 3 are those with very low injection doses. Experimentally it is known that, for example, aCF beam into a SiO2 substrate causes etching for a very shortinitial time period and then results in deposition offluorocarbon films on the substrate in this energy range. Suchexperimental observation is consistent with our MD simulation results given in Fig. 3 although we have not been able
to make quantitative comparison of sputtering yields between
MD simulations and beam experiments yet. Our recent results
based on longer simulation runs indeed indicate the process
turns from etching to deposition as the accumulation ofFig. 3

Acknowledgements
This work was supported by Semiconductor TechnologyAcademic Research Center (STARC).

References[1] S. Hamaguchi, IBM J. Res. Develop. 43, 199 (1999).
[2] S. Hamaguchi, Thin Films 22, 81 (1996).
[3] F.H. Stillinger and T.A. Weber, Phys. Rev. B31, 5262
(1985), J. Chem. Phys. 88, 5123, (1988); Phys. Rev. Lett.
62, 2144 (1989).
[4] H. Ohta and S. Hamaguchi, J. Vac. Sci. Technol. A19,2373 (2001).
[5] H. Ohta and S. Hamaguchi, J. Chem. Phys. 115, 6679
(2001).
[6] H. Yamada and S. Hamaguchi, Proceedings of Symposium on Dry Process, Tokyo, 2002 (IEEJ, TechnicalCommittee on Electron Devices, 2002), 180.
[7] H. Yamada and S. Hamaguchi, J. Appl. Phys. 96, 6174
(2004).
[8] H. Yamada and S. Hamaguchi, J. Plasma Fusion Res.
(present issue).
[9] H. Ohta and S. Hamaguchi, J. Plasma Fusion Res. (present issue).
[10] G.C. Abell, Phys. Rev. B31, 6184 (1985).
[11] J. Tersoff, Phys. Rev. B37, 6991 (1988); Phys. Rev. B38,9902 (1988); Phys. Rev. B39, 5566 (1989).
[12] D.W. Brenner, Phys. Rev. B42, 9458 (1990).

Sputtering yields of SiO2 by CFx (x = 1, 2, 3) and F
beams as functions of injection energy, obtained fromMD simulations. The yield generally depends on the injection dose as well. The data presented here are yields
when injection doses are very low. Van der Waals interactions are not taken into account in these simulations.

483