# SLKMC-II study of self-diffusion of small Ni clusters on Ni (111) surface

Syed Islamuddin Shah<EMAIL_ADDRESS>Giridhar Nandipati
<EMAIL_ADDRESS>Abdelkader Kara<EMAIL_ADDRESS>Talat S.
Rahman<EMAIL_ADDRESS>Department of Physics, University of Central
Florida, Orlando, FL 32816

###### Abstract

We studied self-diffusion of small 2D Ni islands (consisting of up to $10$
atoms) on Ni (111) surface using a self-learning kinetic Monte Carlo (SLKMC-
II) method with an improved pattern-recognition scheme that allows inclusion
of both fcc and hcp sites in the simulations. In an SLKMC simulation, a
database holds information about the local neighborhood of an atom and
associated processes that is accumulated on-the-fly as the simulation
proceeds. In this study, these diffusion processes were identified using the
drag method, and their activation barriers calculated using a semi-empirical
interaction potential based on the embedded-atom method. Although a variety of
concerted, multi-atom and single-atom processes were automatically revealed in
our simulations, we found that these small islands diffuse primarily via
concerted diffusion processes. We report diffusion coefficients for each
island size at various tepmratures, the effective energy barrier for islands
of each size and the processes most responsible for diffusion of islands of
various sizes, including concerted and multi-atom processes that are not
accessible under SLKMC-I or in short time-scale MD simulations.

###### pacs:

68.35.Fx, 68.43.Jk,81.15.Aa,68.37.-d

##  I Introduction

Surface diffusion is of interest not only because it is so different from
diffusion in bulk solids Antczak and Ehrlich (2010) but because diffusion of
adatoms on metal surfaces, individually or as a group via multi-atom or
concerted diffusion processes plays an essential role in a wide variety of
such surface phenomena as heterogenous catalysis, epitaxial crystal growth,
surface reconstruction, phase transitions, segregation, and sintering.Zangwill
(1988) A precise knowledge of diffusion mechanisms is essential for
understanding and control of these phenomena.Kaxiras (1996) Adatoms can
diffuse on a substrate in a variety of ways, and competition between various
types of diffusion processes (due to the differences in their rates)
determines the shapes of the islands formed and (on macroscopic times scales),
the morphological evolution of thin films. Hence a great deal of effort has
been devoted to investigation of self-diffusion of adatom islands on metal
surfaces, initially using field ion microscopy (FIM),Bassett (1976); TSong and
Casanova (1980); Wang and Ehrlich (1990, 1992); Wang  _et al._ (1998); Kellogg
(1993) and more recently scanning tunneling microscopy (STM).Wen  _et al._
(1994); Pai  _et al._ (1997); Giesen and Ibach (2003); Giesen  _et al._
(1998); Fern  _et al._ (2000); der Vegt  _et al._ (1995); Busse  _et al._
(2003); Muller  _et al._ (2005); Liu and Adams (1992) Because of inherent
differences in the microscopic processes responsible for island diffusion on
different metal surfaces, this is still an on-going research problem. Both
experimental and theoretical studies for various systems have succeeded in
finding the activation barriers and prefactors for a single-adatom diffusion
processes Tung and Graham (1980); Fu and Tsong (2000); Flahive and Graham
(1980); Liu  _et al._ (1991); Rice  _et al._ (1992); Stoltze (1994); Li and
DePristo (1996); Mortensen  _et al._ (1996); Chang  _et al._ (2000); Haug and
Jenkins (2000); Kurpick (2001); Bulou and Massobrio (2005); Kim  _et al._
(2006). Ref Antczak and Ehrlich, 2010 provides a good survey of those efforts.
To the best of our knowledge, however, there has so far been no systematic
experimental or theoretical effort to identify the diffusion mechanisms
responsible for diffusion of small 2D Ni islands on Ni(111) and to calculate
their activation barriers. In this article we report our results of doing so
for such islands, ranging in size from 1 to 10 atoms.

Arrangement of atoms in the substrate of an fcc(111) surface results in two
types of three-fold hollow sites for an adatom: the regular fcc site (with no
atom beneath it in the second layer), and an hcp site (with an atom beneath it
in the second layer). Occupancy of adatoms at fcc sites maintains the crystal
stacking order (ABC stacking) of fcc structure, while occupancy of hcp sites
leads to a stacking fault. Depending on its relative occupation energy, which
is material dependent, an adatom can occupy one or the other of these sites.
Which site is preferred on the fcc(111) surface affects the way diffusion and
hence growth proceeds. It is therefore important to understand whether the
diffusion proceeds via movement of atoms from fcc-to-fcc or hcp-to-hcp or fcc-
to-hcp hcp-to-fcc sites. It has been observed experimetally that for smaller
clusters mixed occupancy Repp  _et al._ (2003) of fcc & hcp sites is possible.

A host of studies has been devoted to problems of self-diffusion and diffusion
mechanisms on metal fcc(111) surfaces, almost exclusively, however, with
either a preconceived set of processes or merely approximate activation
barriers. It is nevertheless crucial to discover the full range of processes
at work and to accurately establish the activation barrier of each. It is also
well known that the fcc(111) surface, being atomically flat, has the least
corrugated potential energy surface of any fcc surface, resulting in low
diffusion barriers even for clusters to diffuse as a whole. Consequently,
studies of diffusion processes on fcc(111) surfaces is a challenging problem
for both experiment and simulation even to this day. For a monomer and smaller
islands like dimer, trimer and up to certain extent, tetramer, all possible
diffusion processes may be guessed. But as islands further increase in size,
it becomes more difficult to enumerate all possible diffusion processes a
priori. An alternative is to resort to molecular dynamics (MD) simulation. But
because diffusion processes are rare events, an MD simulation cannot capture
every microscopic process possible, as most of the computational time is spent
in simulating atomic vibration of atoms. Instead, to do a systematic study of
small Ni island diffusion on Ni (111) surface we resorted to an on-lattice
self-learning kinetic Monte Carlo (SLKMC-II) method, which enables us to study
longer time-scales than are feasible with MD yet to find all the relevant
atomic processes and their activation barriers on-the-fly, as KMC methods
limited to a priori set of processes cannot do. Moreover, whereas previous
studies have used an on-lattice SLKMC method, Trushin  _et al._ (2005); Karim
_et al._ (2006); Nandipati  _et al._ (2009, 2011) in which adatoms were
restricted to fcc occupancy, in the present study both fcc and hcp occupancies
are allowed, and are detectible by our recently developed improved pattern-
recognition schemeShah  _et al._ (2012).

The remainder of the paper is organized as follows. In Section II we discuss
the details of our SLKMC-II simulations, with particular attention to the way
we find diffusion processes and calculate their activation barriers. In
Section III we present details of concerted, important multi-atom and single-
atom diffusion processes responsible for the diffusion of Ni islands as a
function of island size. In Section IV we present a quantitative analysis of
diffusion coefficients at various temperatures and of effective energy
barriers as a function of island size. In Section V we present our
conclusions.

##  II Simulation Details

To study Ni island diffusion on fcc Ni(111) surface, we carried out SLKMC
simulations using the pattern-recognition scheme we developed recently Shah
_et al._ (2012) that includes both fcc and hcp sites in the identification of
an atom’s neighborhood. Various types of diffusion processes are possible, and
their activation barrier depends on the atom’s local neighborhood. Whenever a
new neighborhood around an atom is identified, a saddle-point search is
carried out to find all the possible atomic processes and calculate their
activation barriers – provided that it has at least one similar empty site in
the second ringShah  _et al._ (2012), since when an atom occupies an fcc (or
alternatively an hcp) site, the nearest neighbor (NN) hcp (or,
correspondingly, fcc) sites cannot be occupied. In our simulations we used a
system size of 16x16x5 with the bottom 2 layers fixed, and carried out saddle-
point searches using the drag method. In this method a central or active atom
is dragged in small steps towards a probable final position. If the central
atom is on an fcc (hcp) site, then it is dragged towards a NN vacant fcc (hcp)
site in the second ring. Since atoms are allowed to occupy either hcp or fcc
sites, an atom being dragged from an fcc (hcp) site to a neighboring similar
site is allowed to relax to an intermediate hcp (fcc) site in between the two
fcc (hcp) sites. In other words, processes are possible in which atoms in an
island may occupy fcc, hcp or both fcc & hcp sites simultaneously. In the drag
method, the atom being dragged is always constrained in the direction of the
reaction coordinate but allowed to relax along its other degrees of freedom
(those perpendicular to the reaction coordinate), while all the other atoms in
the system are allowed to relax in all degrees of freedom. Once the transition
state is found, the entire system is completely relaxed to find the final
state of the process. The activation barrier of the process is the difference
between the energies of the transition and initial states. We have verified
the activation barriers of some of the key processes found by the drag method
using the (more accurate but computationaly expensive) nudged-elastic band
(NEB) methodJónsson  _et al._ (1998), and found no significant difference. For
inter-atomic interactions, we used an interaction potential based on the
embedded-atom method (EAM) as developed by Foiles et al.Foiles  _et al._
(1986). In all our SLKMC simulations we used the same pre-exponential factor
of $10^{12}$s-1, which has been demonstrated to be a good assumption for such
systems as the one under examination here.Yildirim  _et al._ (2007, 2005)

For the small islands under study here (1-10 atoms), we found that when an
atom is dragged rest of the atoms in the island usually follow. For very small
islands (1-4 atoms), $\it all$ of the processes identified by the drag method
were concerted-diffusion processes. As island size increases we found single-
atom and multi-atom processes as well. For islands of size 5-6, even single-
atom detachment processes are identified and stored in the database (even
though they are not allowed in our simulations). To account for all types of
processes associated with both compact and non-compact shapes – especially
concerted processes and multi-atom processes – we used $10$ rings to identify
the neighborhood around an active atom in our SLKMC simulations. Using $10$
rings corresponds to including fifth nearest-neighbor interactions. To make
sure we identified all the single-atom processes, we also carried out saddle-
point searches with all of the atoms fixed except the atom being dragged.
Although there is no infallible method for discovering all possible processes,
we did exhaust the search for possible processes identifiable using the drag
method.

In order to save computational time, we first carried out SLKMC simulations at
700K for each island size, and used the database thus generated to carry out
our simulations for the same size at lower tempratures (300, 400, 500 and 600K
). The rationale for this approach is that an island goes through many more
shapes at higher tempratures: when a simulation is carried out at a lower
temprature starting out with a database generated at a higher temprature, it
only rarely finds an unknown configuration. It is not possible, however, to
economize on computational time by using, for the smaller islands under study
here, a database generated for (say) the larger among them, because the types
of processes possible (along with their respective barriers) are dependent on
an island’s particular size.

##  III Results

As mentioned above, all of the processes for a given island are identified and
their activation barriers calculated, and stored in a database on-the-fly. We
discuss in this section, however, only key processes of the various general
types (concerted, multi-atom and single-atom).

Fig. 1(a) is a sketch of the fcc(111) surface with its adsorption sites marked
as fcc and hcp. Determining whether an adatom is on an fcc or on an hcp site
on this surface requires knowledge of at least $2$ substrate layers below the
adatom layer. In all our figures we show only the adatom layer and the layer
below (the top substrate layer) with the convention that the center of an
upward-pointing triangle (along the y-axis) formed by the (top layer)
substrate atoms is an fcc site, while the center of a downward-pointing
triangle pointing triangle is an hcp site. An island on an fcc (111) surface
can be on fcc sites or on hcp sites or a combination of both sites (some atoms
of the island sitting on fcc sites and the rest on hcp sites). Depending on
the type of material either the fcc or the hcp site will be energetically
favorable. As we shall see for each island size under study here, the fcc site
for Ni(111) is always at least slightly more favorable than the hcp site.

Figure 1: (a) fcc and hcp sites on an fcc(111) surface, with corresponding
directions for concerted diffusion processes; (b) A-type and B-type step edges
(here, for an all-hcp island) for the same surface.

A compact adatom island on an fcc(111) surface can move in the three
directions shown in Fig.1(a). Note that the numbering scheme for the
directions open to an atom on an fcc site is inverse to that for those open to
an atom on an hcp site (see Fig. 1(a)). We follow the enumeration convention
for directions distinguished in Fig. 1(a) throughout the article in tabulating
activation barriers for concerted processes for islands of various sizes and
shapes. Concerted processes involve all atoms moving together from all-fcc
sites to all-hcp sites or vice-versa. In a concerted diffusion process a
cluster can either translate in one of the three directions shown in Fig. 1
(concerted translation) or rotate around an axis (around the center of mass),
either clockwise or anti-clockwise (concerted rotation). Since concerted
rotational processes do not produce any displacement in the center of mass of
an island, they do not contribute to island diffusion. Depending on the size
of the island and its shape, activation barriers for the processes in these
three directions can be different.

Activation barriers for single-atom processes, however, depend on the type of
step-edge along which atom diffuses. Fig. 1(b) shows, using the example of a
6-atom hcp island, how an A-type step-edge $\\--$ a (100) micro-step differs
from a B-type step-edge $\\--$ a (111) micro-step. We discuss important
single-atom diffusion processes systematically and in detail in Sub-section
III.11.

As island size increases not only does the frequency of single-atom processes
increase but the frequency of multi-atom processes does so as well. All multi-
atom mechanisms involve shearing. A special case is reptation mechanismChirita
_et al._ (2000, 1999), a two-step shearing process that moves the cluster from
all-fcc to all-hcp sites or the reverse: first, part of the island moves from
fcc to hcp sites; then the rest of the island moves from fcc to hcp. Hence at
the intermediate stage, the island has mixed fcc-hcp occupancy. In case of Ni-
island diffusion, reptation processes occur only when the shape of the island
becomes non-compact. We will discuss reptation in detail when we take up
islands of size 8-10.

###  III.1 Monomer

As mentioned earlier, much work has been done to determine activation barrier
for Ni monomer diffusion on Ni (111) surfaceAntczak and Ehrlich (2010). A
monomer on fcc (111) surface can adsorb either on an fcc or an hcp site. We
find that adsorption of an adatom on an fcc site is slightly favored over than
on an hcp site by $0.002$ eV – in good agreement with the value reported in
Ref Liu  _et al._ , 1991. Diffusion of a monomer occurs through hopping
between fcc sites via an intermediate hcp site. We find the activation barrier
for a monomer’s hopping from an fcc site to a neighboring hcp site to be
$0.059$ eV while that for the reverse process is $0.057$ eV. The effective
energy barrier for monomer is found to be 0.057 eV, which is consistent with
the result reported by Liu et al.Liu and Adams (1992) of $0.056$ eV.

###  III.2 Dimer

On any fcc(111) surface a dimer (of the same species) can have three possible
arrangements: both atoms on fcc sites (an FF-dimer, Fig. 2(b)), both on hcp
sites (an HH-dimer, Fig. 2(a)) or one atom on an fcc and the other on an hcp
site (an FH-dimer, Fig. 2(c)). We find that the FF-dimer is energetically more
favorable than the HH-dimer by $0.005$eV and the FH-dimer the least favorable
by $0.011$ eV. We find that both FF and HH dimers diffuse via concerted as
well as single-atom processes, whereas the FH-dimer diffuses via single-atom
processes only.

Figure 2: Possible configurations for a dimer, with activation barriers for concerted diffusion processes. (a) FF dimer (both atoms on fcc sites); (b) HH dimer (both atoms on hcp sites); (c) FH dimer (one atom on an fcc and the other on an hcp site); (d) FF dimer in concerted clockwise rotation and (e) HH dimer in concerted counter-clockwise rotation. Table 1: Activation barriers (in eV) of concerted processes for dimer diffusion. Direction | fcc | hcp
---|---|---
1 | 0.071 | 0.066
2 | 0.148 | 0.143
3 | 0.148 | 0.143

In concerted diffusion processes, both atoms in a FF (HH) dimer move from fcc
(hcp) to the nearest hcp (fcc) sites as shown in Fig. 2 (a)&(b) and Fig. 2
(d)&(e), thereby converting an FF (HH) dimer into an HH (FF) dimer. In the
case of an FF (HH) dimer, the activation barrier for concerted translational
(Fig. 2 (a)&(b)) is $0.148$ eV ($0.143$ eV) while that for concerted rotation
(Fig. 2 (d)&(e)) is $0.038$ eV. In concerted dimer rotation, the activation
barriers for both clockwise and anti-clockwise directions are the same, as
they are symmetric to each other. Activation barriers for translational
concerted diffusion processes in all three directions (see Fig. 1) both for FF
and HH dimers are reported in Table. 1. Our results for concerted processes
are $0.028$ eV higher than the corresponding activation barriers for a dimer
reported in Ref. Liu and Adams, 1992. (This difference $\\--$ as with those in
what follows $\\--$ may be due to the different inter-atomic potential
employed in their study and ours.)

Single-atom processes transform both FF and HH dimers into an FH-dimer. In
this case one of the fcc atoms in an FF-dimer or an hcp atom in an HH-dimer
moves to a nearest-neighbor hcp or an fcc site respectively, as shown in Fig.
2 (a)&(b) with the double-headed arrow. The activation barriers are $0.034$ eV
and $0.035$ eV for hcp and fcc dimer, respectively. In the case of an FH-
dimer, two types of single-atom diffusion processes are possible, as shown in
Fig. 2: an fcc atom moves to the nearest hcp site in the direction of the open
arrowhead, forming an HH-dimer, or an hcp atom moves in the direction of the
solid arrowhead to the nearest fcc site, forming an FF-dimer. The activation
barriers for these processes are $0.028$ eV and $0.024$ eV, respectively.

###  III.3 Trimer

Depending on where a third atom is attached to the dimers shown in Fig. 2(a &
b), there are four possible arrangements of atoms in a compact trimer: two
types of fcc timers – one centered around an hcp site (F3H), the other
centered around a top site (F3T) (see Figs. 3(a) & (d)) $\\--$ and two types
of hcp trimers $\\--$ one centered around an fcc site (H3F), the other
centered around a top site (H3T) (Figs. 3(c) & (b)). Although all four trimers
have the same shape, their local environment is different, so that their
adsorption energies are distinct, as are the activation barriers for their
possible diffusion processes. F3T timer is the most energetically favorable:
F3H, H3T and H3F are less energetically favorable by $0.006$, $0.007$ and
$0.0013$ eV, respectively. It should also be noted that although trimers can
take on non-compact shapes, the configurations depicted in Fig. 3 are the most
frequently observed in our trimer simulations.

Figure 3: Possible arrangements of atoms in a trimer, with possible concerted
diffusion processes and their activation barriers. (a)$\\--$(d) Concerted
translation: (a) F3H-all atoms on fcc sites centered around an hcp site; (b)
H3T-all atoms on hcp sites centered around a top site; (c) H3F-all atoms on
hcp sites centered around an fcc site; (d) F3T-all atoms on fcc sites centered
around a top site. (e) & (f) Concerted rotation: F3T and H3T respectively.

In the case of F3T and H3T trimers, two types of concerted processes were
observed, a non-diffusive concerted rotation (clockwise and anti-clockwise)
(Fig 3 (e) & (f)) and a diffusive concerted translation (in all three
directions) (Fig 3 (d) & (b)). Concerted rotation processes transform an H3T
timer into a F3T trimer and vice versa. The activation barrier for the
concerted rotation processes for F3T trimer is $0.150$ eV while that for those
of the H3T trimer is $0.144$ eV. Translation transforms an F3T timer into a
H3F timer and vice versa. The activation barrier for the concerted
translations possible for these two trimers are $0.200$ eV and $0.193$ eV for
F3T and H3T, respectively. The activation barrier of $0.200$ eV for
translational motion of F3T trimer is in agreement with the value reported in
Ref Liu and Adams, 1992. For F3H and H3F trimers, only concerted translation
processes are possible; their activation barriers are $0.194$ eV and $0.186$
eV (the value reported for the same process in Ref Liu and Adams, 1992 is
$0.187$ eV), respectively. Fig. 3 (a-b) & (c-d) reveal that these concerted
diffusion processes transform an F3H into an H3T timer and an H3F to an F3T
trimer. Since the shape of these trimers is symmetric (see Fig 3 (a) & (c)),
the activation barriers for their diffusion in all $3$ possible directions are
the same.

Table 2: Activation barriers (in eV) for single-atom diffusion processes for an H3T compact trimer in the directions shown in Fig. 4.  Type | 1 | 2 | 3 | 4
---|---|---|---|---
F3T | 0.439 | 0.858 | 0.858 | 0.439
F3H | 0.432 | 0.875 | 0.875 | 0.432
H3T | 0.436 | 0.856 | 0.856 | 0.436
H3F | 0.429 | 0.872 | 0.872 | 0.429

As for single-atom processes in the case of a trimer: an atom can move in $4$
different directions as shown in Fig. 4, resulting $2$ different types of
single-atom processes: directions $1$ & $4$ correspond to edge-diffusion
processes which open up the trimer; directions $2$ & $3$ correspond to
detachment processes (excluded from the present study, which is confined to
diffusion of single whole islands, in which an island’s integrity [and hence
its size] is maintained). We note that these processes move atoms from fcc
(hcp) to nearest fcc (hcp) site. Activation barriers for processes in these
$4$ directions for different types of trimers are given in Table. 2. Because
these activation barriers are so high relative to those for concerted
processes, single-atom processes were rarely observed in our simulations of
trimer diffusion.

Figure 4: Single-atom processes possible for an H3T trimer. Activation
barriers for the processes in these 4 directions for the 4 posible trimer
configurations are given in Table. 2

We note that as island size increases, possible types of single-atom processes
increases as well (though with the decamer, basically all possible types have
appeared). Accordingly, it is convenient to defer detailed discussion of
single-atom processes until later (Sub-section. III.11)

###  III.4 Tetramer

Adding another atom to any of the trimers shown in Fig. 3 (a - d) results in
the formation of a compact tetramer, diamond-shaped $\\--$ with a long
diagonal (along the line joining farthest atoms) and a short one
(perpendicular to the long one), as shown in Fig. 5. Once again the fcc island
(Fig. 5(a)) is energetically more favorable than the hcp one $\\--$ in this
case, by $0.009$ eV. Three types of translational concerted diffusion
processes are possible for each of the fcc and hcp tetramers, that is, one
along each of the three directions specified in Fig.1. An example of a
cncerted fcc-to-hcp process (along direction 1) for a tetramer is shown in
Fig. 5(a); its activation barrier is 0.213 eV. The reverse process (hcp to
fcc) is shown in Fig. 5(b); its activation barrier is 0.204 eV (the value
reported in Ref. 18 is 0.210 eV). Because the process in direction 3 is
symmetric to that in direction 1, the energy barriers for these processes are
identical, as are those for the reverse processes. The energy barrier along
direction 2 is 0.313 eV from fcc to hcp and 0.304 eV from hcp to fcc. These
values are systematically displayed in Table. 3.

Although the multi-atom processes shown in Fig. 5(c & d) have activation
barriers lower than those of single-atom processes, they could not be found
using the drag method. These processes thus had to be added manually to the
database and their activation barriers obtained using the nudged elastic band
method (NEB). In these multi-atom processes, two atoms move togehter in the
same direction, the result is a shearing mechanism as shown in Fig. 5(c) &
(d). For this shearing process, from fcc to hcp, the activation barrier is
0.285 eV; that for the reverse process from hcp to fcc is 0.276 eV. The drag
method also finds single-atom processes, but because in tetramers (as in
islands of size 3 $\\--$ 7) these have higher activation barriers than those
of concerted processes, they were not observed during the simulations.

Figure 5: Diffusion processes possible for a tetramer: (a) & (b) concerted diffusion along the short diagonal; (c) & (d) shearing processes. Table 3: Activation barriers (eV) for the concerted tetramer translation processes shown in Fig. 5(a)&(b).  Direction | fcc | hcp
---|---|---
1 | 0.213 | 0.204
2 | 0.313 | 0.304
3 | 0.213 | 0.204

###  III.5 Pentamer

Figure 6: Examples of concerted diffusion processes for a compact pentamer,
along with their activation barriers.

The compact shapes of a pentamer can be obtained by attaching an atom to a
diamond-shaped tetramer. Although the geometries of compact pentamer clusters
thus obtained are the same, the island’s diffusion is crucially affected by
where this additional atom is placed: attachment of an atom to an A-type step-
edge of an fcc tetramer results in the long A-type step-edge pentamer shown in
Fig. 6(c); attachment of an atom to a B-type step-edge of the same tetramer
results in the long B-type step-edge pentamer shown in Fig. 6(a)); the
corresponding results of attaching an atom to an hcp tetramer are shown in
Figs. 6(b) and (d), respectively. The most energetically favorable of these is
the fcc pentamer with a long A-type step-edge (Fig. 6(c)); less favorlable by
$0.005$ eV are the fcc pentamer with a long B-type step-edge (Fig. 6(a)), by
$0.011$ eV the hcp pentamer with a long A-type step-edge, and by $0.017$ eV
the hcp pentamer with a long B-type step-edge. That is: as usual, fcc islands
are more stable than hcp ones. And, within each of those types, pentamers with
a long A-type step-edge are more stable than those with a long B-type step-
edge.

Figure 7: The single-atom processes that convert the long A-type step-edge pentamer to the long B-type step-edge pentamer. Table 4: Activation barriers (eV) without parantheses are for concerted-translation processes of pentamers with a long A-type step-edge, as shown in Figs. 6 (c)&(b); barriers in parentheses are for such processes for pentamers with a long B-type step-edge, as shown in Figs. 6 (a)&(d). Directions | fcc | hcp
---|---|---
| A (B) | A (B)
1 | 0.348 (0.301) | 0.342 (0.284)
2 | 0.348 (0.353) | 0.342 (0.337)
3 | 0.295 (0.353) | 0.289 (0.337)

In our simulations we found that compact pentamers diffuse mostly via
concerted diffusion processes, which displace the island as a whole from fcc-
to-hcp or vice-versa. Fig. 6 shows concerted diffusion processes along
direction 1 for the long B-type step-edge pentamer and along direction 3 for
the long A-type step-edge pentamer. Table 4 displays activation barriers for
concerted processes in all 3 directions for both types of pentamer. Fig. 7
(a),(b) & (c) shows the single-atom processes that transform an fcc pentamer
from long A-type (= short B-type) to a short A-type (= long B-type) cluster,
with the activation barrier for each.

###  III.6 Hexamer

Depending on whether a sixth atom is attached to a long A-type or to a long
B-type step-edge pentamer, there are $3$ possible compact shapes for a
hexamer: (1) when an atom is added in such a way as to extend the shorter edge
of either a long A-type or a long B-type step-edge pentamer, the result is one
of the parallelogramic hexamers shown in Fig. 8; (2) when an atom is attached
to the long edge of either type of pentamer, the result is the one of the
irregular hexamers shown in Fig. 9; ($3$) when an atom is added to the shorter
edge of either type of pentamer, the result is one of the triangular hexamers
with all step edges of either the A-type (Figs. 10(b) &(c)) or of the B-type
(Figs. 10(a) &(d)).

Figure 8: Parallelogramic hexamers obtained by extending (by one atom) the shorter edge of either the long A-type or long B-type pentamers: (a) fcc cluster; (b) hcp cluster. The activation barriers indicated are for concerted diffusion in direction $1$. Table 5: Relative stabilities of the hexamers most frequently observed in our simulations. P = Parallelogram; I = Irregular; T = Triangular. Type | Shape | Description | Energy (eV) | Reference
---|---|---|---|---
fcc | P | equal A & B steps | 0 | Fig. 8(a)
fcc | I | edge atom on A step | 0.011 | Fig. 9(b)
fcc | I | edge atom on B step | 0.011 | not shown
hcp | P | equal A & B steps | 0.014 | Fig. 8(b)
fcc | T | all A steps | 0.019 | Fig. 10(c)
hcp | I | edge atom on B step | 0.024 | Fig. 9(b)
hcp | I | edge atom on A step | 0.024 | not shown
fcc | T | all B steps | 0.029 | Fig. 10(a)
hcp | T | all A steps | 0.032 | Fig. 10(b)
hcp | T | all B steps | 0.042 | Fig. 10(d)

Table. 5 shows the order of relative stabilties of the hexamers most
frequently obseved in our simulations. It reveals that, for $\it a$ $\it
given$ $\it shape$, hexamers on fcc sites are more energetically favored than
those on hcp sites, and that among hexamers on fcc sites (as for those on hcp
sites), clusters in which A steps are longer than B steps are more stable.

Figs. 8-10 also show concerted diffusion processes (in direction 1) for these
hexamers, together with the activation barriers for each. Tables. 6 & 7 give
activation barriers for the hexamers shown in Figs. 8 & 9, respectively. Since
triangular hexamers (Fig. 10) are symmetric, their activation barriers for
concerted diffusion are same in all three directions.

Figure 9: Irregular hexamers obtained by attaching an atom to the long edge of a pentamer: (a) fcc cluster; (b) hcp cluster. The activation barriers indicated are for concerted diffusion in direction $1$. Figure 10: Triangular hexamers obtained by adding an atom to the short edge of a pentamer: (a) fcc hexamer with B-type step edges; (b) hcp hexamer with A-type step edges; (c) fcc hexamer with A-type step edges; (d) hcp hexamer with B-type step edges. The activation barriers here are for concerted diffusion in direction $1$. Table 6: Activation barriers (eV) of the concerted translation processes in all three directions for the hexamer shown in Fig.8 Directions | fcc | hcp
---|---|---
1 | 0.374 | 0.360
2 | 0.466 | 0.451
3 | 0.254 | 0.240
Table 7: Activation barriers (eV) of concerted translations processes in all $3$ directions for the hexamers shown in Fig.9 Directions | fcc | hcp
---|---|---
1 | 0.393 | 0.380
2 | 0.397 | 0.383
3 | 0.391 | 0.378

Fig. 11 shows the most frequently observed multi-atom processes for a hexamer
– shearing processes in which a dimer moves along the A-type step-edge of the
cluster from sites of one type to the nearest-neighbor sites of the same type.
Figs. 11 (a) & (b) show this kind of diffusion process for an hcp cluster and
Figs. 11 (c) & (d) for an fcc cluster. Although this dimer shearing process
does not much displace the center of mass of a hexamer; it does have a
striking consequence: it converts a parallelogramic hexamer (Fig. 8) into an
irregular hexamer, (Fig. 9) and vice-versa.

Figure 11: Dimer shearing processes in case of a hexamer along with their
activation barriers: (a) & (b) all-hcp hexamers; (c) & (d) all-fcc clusters.

###  III.7 Heptamer

Figure 12: Concerted diffusion processes and their activation barriers in
direction $1$ for a heptamer.

On an fcc(111) surface, an heptamer has a compact closed-shell structure with
each edge atom having at least three nearest-neighbor bonds, as shown in Fig.
12. Our SLKMC simulations (keep in mind here the range of temperatures to
which they were confined) found that heptamer diffuses exclusively via
concerted diffusion processes, which displace the cluster from fcc-to-hcp and
vice versa; the barriers for which are shown in Fig. 12. That these processes
will predominate can also be concluded from the fact that the effective energy
barrier for heptamer diffusion (cf. Table. 23) is close to the average of the
activation barriers shown in Fig. 12.

Since the compact heptamer has a symmetric shape, activation barriers in all
three directions are the same as those shown in Fig. 12. Again, the fcc island
is more energetically favorable than its hcp counterpart – in this case by
$0.015$ eV.

###  III.8 Octamer

Figure 13: Possible orientations for a compact fcc octamer: (a) with long
A-type step edges; (b) with long B-type step edges.

Compact octamers have two distinct orientations, one with two long A-type
step-edges, the other with two long B-type step-edges, as shown in Fig. 13 (a)
& (b) for an fcc octamer. Octamers with long A-type step-edges (Fig. 13 (a))
can be obtained by attaching an atom to any B-type step-edge of a compact
heptamer, while a compact octamer with long B-type step-edges results from
attaching one to any A-type step-edge. Again: the fcc islands are more
energetically favorable than the hcp ones, and within each type, islands with
long A-type step edges are more stable than those with long B-type step-edges.

A compact octamer diffuses via concerted diffusion processes, as shown in Fig
14. The activation barrier of a concerted diffusion process depends on whether
the octamer has long A- or long B-type step-edges. As Fig. 14 shows, a
concerted diffusion process converts a long A-type step-edge fcc octamer into
a long B-type step-edge hcp cluster, and vice-versa. Table. 8 reports the
activation barriers for concerted diffusion processes in all $3$ directions
for both orientations for an fcc as well as an hcp (see Fig. 13) octamer.
Concerted diffusion processes in directions $2$ and $3$ are the most
frequently observed processes in octamer diffusion.

Figure 14: Examples of concerted diffusion processes in direction 2 for an octamer: (a) compact fcc octamer with long A-type step-edges; (b) compact hcp octamer with long B-type step-edges. Table 8: Activation barriers (in eV) of the twelve concerted diffusion processes for compact octamers. Directions | fcc | hcp
---|---|---
| A (B) | A (B)
1 | 0.589 (0.585) | 0.567 (0.571)
2 | 0.491 (0.484) | 0.468 (0.469)
3 | 0.491 (0.484) | 0.468 (0.469)

Although an octamer diffuses primarily via concerted processes we found in our
simulations that both multi-atom and single-atom processes are also relatively
common. As mentioned before, we defer comprehensive discussion of single-atom
processes to section III.11. Here we discuss multi-atom processes particular
to octamers.

Figure 15: Example of shearing diffusion processes within an octamer along
with their activation barriers.

Multi-atom processes involving shearing and reptation are shown in Figs. 15 &
16, respectively. In shearing processes, $\it part$ of the island (more than
one atom) moves from fcc to the nearest fcc sites, if all the island is
initially on fcc sites – or from hcp to hcp sites, if all of it initially sits
on hcp sites. Figs. 15 (a) & (b) show trimer shearing processes within an
octamer, along with their activation barriers, while Figs. 15 (c) & (d) show
dimer shearing processes, with their activation barriers.

Reptation is a 2-step diffusion process. In the case of an fcc island, the
entire island diffuses from fcc to nearest-neighbor hcp sites in two steps. In
the first, part of the island moves from fcc to nearest-neighbor hcp sites,
leaving part of the island on fcc sites and part on hcp sites. In the next,
the remainder of the island initially on fcc sites moves to hcp sites. Figs.
16 (a)-(d) shows various steps (subprocesses) of a reptation process, with
their activation barriers.

Figure 16: Successive sub-processes (or steps) involved in an octamer
reptation diffusion mechanism.

###  III.9 Nonamer

Figure 17: Concerted diffusion processes and their activation barriers for
nonamers. Figure 18: Dimer shearing processes and their activation barriers
for compact nonamers.

For a nonamer, we observed all types of diffusion processes – single-atom,
multi-atom and concerted. The most frequently observed were two types of
single-atom mechanisms: edge-diffusion processes along an A- or a B-type step-
edge and corner rounding (Fig. 21). The nonamer is the smallest island for
which concerted processes are not the most frequently picked (the concerted
processes for the most frequently observed nonamer configurations [compact or
nearly so] are shown in Fig. 17). Even so, concerted processes contribute the
most to $\it island$ diffusion: that is, the displacement they produce in the
nonamer’s center of mass is far greater than that produced by single-atom
processes, despite the far greater frequency of the latter. This is reflected
in the fact that the effective energy barrier for nonamer (cf. Table 12) is
much closer to the average activation barrier for concerted processes (cf.
Table. 9) than for that of single-atom processes (cf. Table 9). The fact that
the effective activation barrier is slightly higher than the average energy
barrier for concerted processes is due mainly to the contribution of kink
processes, which do contribute somewhat to island diffusion.

Table 9: Activation barriers (eV) of concerted translations processes in all
$3$ directions for the nonamers shown in Figs.17(a–d).

Directions | fcc | hcp
---|---|---
| (a) | (b)
1 | 0.520 | 0.499
2 | 0.626 | 0.605
3 | 0.605 | 0.583

fcc | hcp
---|---
(c) | (d)
0.486 | 0.465
0.486 | 0.465
0.693 | 0.672

The most frequently observed multi-atom processes are the four forms of dimer
shearing along an A-type step-edge shown in Fig. 18(c)-(f), as have been
discussed above for island of sizes 6 & 8\. The activation barriers for these
dimer shearing processes are lower than those for single-atom diffusion
processes along an edge and also for some corner rounding processes. Reptation
processes also show up, but only when the nonamer is non-compact (we do not
illustrate these here) Shah  _et al._ (2012).

###  III.10 Decamer

Even in the case of a decamer, we have observed single-atom, multi-atom and
concerted diffusion processes. Single-atom diffusion processes are the most
frequently observed. The most frequently observed compact shape of decamer
during our simulations is that shape shown in Fig. 19, which has the same
number of A- and B-type step edges. As usual, an fcc cluster is energetically
more favorable than an hcp cluster.

Table 10: Activation barriers (eV) of concerted diffusion processes in all $3$ directions of decamer as shown in Fig.19 (a) & (b). Directions | fcc | hcp
---|---|---
1 | 0.661 | 0.638
2 | 0.700 | 0.677
3 | 0.700 | 0.677

For the shape shown in Fig. 19, the most frequently observed concerted
diffusion processes are those shown in the same figure, along with their
activation barriers reported in Table 10. It can be seen from the Table 23,
that the effective energy for decamer diffusion is close to that of the
average energy barrier of these concerted processes. That is why decamer
diffusion is dominated by concerted processes.

Figure 19: Frequent concerted diffusion processes and their activation
barriers for compact decamers. Figure 20: Diffusion steps (sub-processes) in
decamer reptation. Note that the decamers in (b) and (c) are identical, though
the arrows indicate different processes.

As with the nonamer, a decamer also undergoes multi-atom processes (shearing
and reptation). Of these, the most frequently observed is dimer shearing along
an A-type step-edge, similar to what has been discussed for clusters of size
6, 8 and 9. Fig. 20 shows the sub-processes in a reptation process, along with
the activation barrier of each.

###  III.11 Single-atom Processes

In this section we provide detail about single-atom processes: edge-diffusion,
corner rounding, kink attachment, and kink detachment, as shown in Fig. 21 for
an hcp island. Their corresponding activation barriers and those for their fcc
analogues are given in Table 11. In each single-atom process, an atom on an
fcc site moves to a nearest-neighbor fcc site, while an atom on hcp site moves
to a nearest-neighbor hcp site. The activation barriers for single-atom
processes depend not only on whether the atom is part of an fcc sland or an
hcp island, but also on whether the diffusing atom is on an A-type or a B-type
step-edge.

In classifying single-atom processes in Table 11 we have used the notation
X${}_{n_{i}}$U $\rightarrow$ Y${}_{n_{f}}$V, where where X or Y = A (for an
A-type step-edge) or B (for a B-type step-edge) or K (for kink) or C (for
corner) or M (for monomer); ni = the number of nearest-neighbors of the
diffusing atom before the process; nf = the number of that atom’s nearest
neighbors after the process. U or V = A or B (for corner or kink processes) or
null (for all other other process types).

For example, process 1, B2 $\rightarrow$ B2, is a single-atom B-step edge
process in which the diffusing atom has 2 nearest-neighbors before and after
the process. Process 3, C1B $\rightarrow$ B2, is a corner rounding process
towards a B-step, the diffusing atom starting on the corner of a B-step with
one nearest-neighbor and ending up on the B-step with two nearest-neighbors.
In process 10, C2A $\rightarrow$ C1B, the diffusing atom begins on the corner
of an A-step having two nearest-neighbors and ends up on the corner of a
B-step with only one nearest-neighbor.

Figure 21: Single-atom processes for an hcp island (Though analogous processes occur for an fcc island, we do not illustrate them here). The index numbers designate the processes described in Table 11, which gives the activation barriers for each. Table 11: Activation barriers (eV) of single-atom processes for both fcc and hcp islands. The index numbers refer to the types of processes illustrated in Fig. 21. See text for explaination of the notation used to classify the process types. Index no. | Process type | fcc | hcp
---|---|---|---
1 | B2 $\rightarrow$ B2 | 0.454 | 0.448
2 | B2 $\rightarrow$ M | 0.821 | 0.815
3 | C1B $\rightarrow$ B2 | 0.177 | 0.173
4 | C1B $\rightarrow$ M | 0.458 | 0.455
5 | C1B $\rightarrow$ A2 | 0.040 | 0.038
6 | C1B $\rightarrow$ C1B | 0.540 | 0.585
7 | C1B $\rightarrow$ M | 0.811 | 0.809
8 | A2 $\rightarrow$ M | 0.795 | 0.794
9 | A2 $\rightarrow$ A2 | 0.326 | 0.307
10 | C2A $\rightarrow$ C1B | 0.399 | 0.397
11 | C2A $\rightarrow$ M | 0.787 | 0.785
12 | B2 $\rightarrow$ K3B | 0.415 | 0.298
13 | K3A $\rightarrow$ A2 | 0.601 | 0.701
14 | A2 $\rightarrow$ K3A | 0.302 | 0.389
15 | K3B $\rightarrow$ B2 | 0.729 | 0.597
16 | K3B $\rightarrow$ C1B | 0.731 | 0.787
17 | K3B $\rightarrow$ M | 1.138 | 1.150
18 | K3B $\rightarrow$ C1B | 0.820 | 0.759

##  IV Diffusion coefficients and effective energy barriers

We start our SLKMC simulations with an empty database. Every time a new
configuration (or neighborhood) is turned up, SLKMC-II finds on the fly all
possible processes using the drag method, calculates their activation barriers
and stores them in the database as the simulation proceeds. Calculation of
energetics occurs at each KMC step during initial stages of the simulation
when the database is empty or nearly so, and ever less frequently later on.
Recall that the types of processes and their activation barriers are dependent
on island size, each one of which requires a separate database that cannot be
derived from that for islands of other sizes.

Figure 22: Arrhenious plots for 1$\\--$10 atom islands.

We carried out 107 KMC steps for each island size at temperatures 300K, 400K,
500K, 600K and 700K. We calculated the diffusion coefficient of an island of a
given size using Einstein Equation:Einstein (1905) $D=\lim_{t\to\infty}\langle
R_{CM}(t)-R_{CM}(0)]^{2}\rangle/2dt$, where $R_{CM}(t)$ is the position of the
center of mass of the island at time $t$, and $d$ is the dimensionality of the
system, which in our case is $2$. Diffusion coefficients thus obtained for
island sizes $1-10$ at various temperatures are summarized in Table 12. At
$300$ K, diffusion coefficients range from $1.63\times 10^{11}$ Å2/s for a
monomer to $8.66\times 10^{01}$ Å2/s for a decamer. Effective energy barriers
for islands are extracted from their respective Arrhenius plots (Fig. 22) and
also summarized in Table 12. Fig. 23 plots effective energy barrier as a
function of island size. It can be seen that the effective energy barrier
increases almost linearly with island size.

Figure 23: Effective energy barriers of 1$\\--$10 atom islands as a function of island size. Table 12: Diffusion coefficients ((Å2/s) at various temperatures and effective energy barriers for Ni islands. Island Size | 300K | 400K | 500K | 600K | 700K |  $E_{eff}$ (eV)
---|---|---|---|---|---|---
1 | $1.63\times 10^{11}$ | $2.85\times 10^{11}$ | $3.99\times 10^{11}$ | $5.00\times 10^{11}$ | $5.87\times 10^{11}$ | 0.058
2 | $8.14\times 10^{09}$ | $2.88\times 10^{10}$ | $6.54\times 10^{10}$ | $1.13\times 10^{11}$ | $1.64\times 10^{11}$ | 0.136
3 | $1.09\times 10^{09}$ | $7.42\times 10^{09}$ | $2.30\times 10^{10}$ | $5.00\times 10^{10}$ | $8.23\times 10^{10}$ | 0.196
4 | $7.64\times 10^{08}$ | $5.83\times 10^{09}$ | $2.02\times 10^{10}$ | $4.91\times 10^{10}$ | $9.93\times 10^{10}$ | 0.217
5 | $2.90\times 10^{06}$ | $8.42\times 10^{07}$ | $6.59\times 10^{08}$ | $2.55\times 10^{09}$ | $7.36\times 10^{09}$ | 0.353
6 | $1.08\times 10^{06}$ | $7.39\times 10^{07}$ | $5.66\times 10^{08}$ | $2.69\times 10^{09}$ | $8.12\times 10^{09}$ | 0.400
7 | $6.24\times 10^{04}$ | $6.40\times 10^{06}$ | $1.01\times 10^{08}$ | $6.37\times 10^{08}$ | $2.37\times 10^{09}$ | 0.477
8 | $1.38\times 10^{04}$ | $1.50\times 10^{06}$ | $2.31\times 10^{07}$ | $1.57\times 10^{08}$ | $5.91\times 10^{08}$ | 0.482
9 | $5.48\times 10^{02}$ | $8.77\times 10^{04}$ | $1.89\times 10^{06}$ | $1.35\times 10^{07}$ | $6.23\times 10^{07}$ | 0.525
10 | $8.66\times 10^{01}$ | $2.25\times 10^{04}$ | $7.20\times 10^{05}$ | $7.97\times 10^{06}$ | $4.95\times 10^{07}$ | 0.627

##  V Further Discussion And Conclusions

To summarize: we have performed a systematic study of the diffusion of small
Ni islands (1-10 atoms) on Ni(111), using a self-learning KMC method with a
newly-developed pattern recognition scheme (SLKMC-II) in which the system is
allowed to evolve through mechanisms of its choice on the basis of a self-
generated database of single-atom, multiple-atom and concerted diffusion
processes (each with its particular activation barrier) involving fcc-fcc,
fcc-hcp and hcp-hcp jumps. We find that concerted diffusion processes
contribute the most to the displacement of the center of mass (i.e., to island
diffusion), while single-atom processes contribute the least. As for multi-
atom processes (reptation or shearing): while these produce more displacement
than the latter, they are hardly ever selected. Though the energy barriers for
reptation processes are small compared to those for concerted diffusion
processes, reptation occurs only when an island is transformed into a non-
compact shape, as happens only rarely in the temprature range to which our
study is confined. In contrast, though shearing occurs with close-to-compact
shapes (which appear more frequently than do non-compact shapes, but with
islands of certain sizes), the barriers for these processes are higher than
those for reptation. Finally, although for all island sizes, island diffusion
is primarily dominated by concerted diffusion processes, the frequency of
occurrence of both single-atom and multi-atom processes does increase with
increase in island size, owing to increase in the activation barrier for
concerted diffusion processes with island size.

Figure 24: Distribution of single-atom, multi-atom, concerted and total
processes for 1-10 atom islands accumulated in the database during SLKMC
simulations. Inset shows the log-linear plot for up to the 6-atom island.

Fig. 24 shows the number of each type of diffusion process (single-, multi-
atom and concerted) collected during our SLKMC-II simulations for each island
size (1-10), together with the total number of processes of all types for each
island size. (For the sake of clarity, the insert shows a log-linear plot of
these quantities for island sizes 1-6). It can be seen that the number of
processes accumulated increases with island size, and significantly so beyond
the tetramer. It can also be seen from Fig. 24 that the overall increase in
number of processes with island size is constituted predominantly by
significant increases in single-atom processes and (to a lesser degree) multi-
atom processes. Meanwhile the number of concerted processes accumulated in the
database increases at a much slower pace with island size.

This significant increase in single-atom processes is mainly due to the use of
$10$ rings to identify the neighborhood of an atom. Use of $10$ rings
corresponds to inclusion of $5$ nearest-neighbor fcc-fcc or hcp-hcp
interactions. Elsewhere we show that $6$ rings (which corresponds to $3$
nearest-neighbor interactions) offer a range of interaction sufficient for
accurately calculating the activation barriers for single-atom processes.Shah
_et al._ (shed) But it is essential to include the long-range interaction (and
hence 10 rings) if one aims to accurately take into account multi-atom and
concerted processes, the latter of which predominate in small-island
diffusion. This significant increase in the number of processes with island
size also justifies resorting to an automatic way of finding all the possible
processes during simulations instead of using a fixed (and thus necessarily
preconceived) list of events.

As mentioned earlier, with increasing island size not only does the number of
accumulated single-atom processes and multi-atom processes increase but also
their frequency of occurrence. Still island diffusion is primarily due to
concerted diffusion processes, since it is these that produce largest
displacement of center of mass. This can be easily observed by comparing, for
each island size, the effective diffusion barriers given in Table 12 with the
activation barriers in the tables given in Sect. III for concerted diffusion
processes: effective diffusion barriers more or less closely follow activation
barriers for concerted diffusion processes – except for the 9-atom island, in
which the contribution of single-atom processes to the island’s diffusion is
significantly larger than for other island sizes.

In conclusion, Ni small-island diffusion on Ni (111) is primarily due to to
concerted diffusion processes, even though the frequency of their occurrence
decreases with increase in island size.

###### Acknowledgements.

We would like to acknowledge computational resources provided by University of
Central Florida. We also thank Lyman Baker for critical reading of the
manuscript and Oleg Trushin for his suggestions during initial stages of this
work.

## References

  * Antczak and Ehrlich (2010) G. Antczak and G. Ehrlich,  _Surface Diffusion : metals, metal atoms and clusters_ (Cambridge University Press, 2010).
  * Zangwill (1988) A. Zangwill,  _Physics at Surfaces_ (Cambridge University Press, 1988).
  * Kaxiras (1996) E. Kaxiras, Comput. Mater. Sci 6, 158 (1996).
  * Bassett (1976) D. W. Bassett, J. Phys. C 9, 2491 (1976).
  * TSong and Casanova (1980) T. T. TSong and R. Casanova, Phys. Rev. B 22, 4632 (1980).
  * Wang and Ehrlich (1990) S. C. Wang and G. Ehrlich, Surf. Sci. 239, 301 (1990).
  * Wang and Ehrlich (1992) S. C. Wang and G. Ehrlich, Phys. Rev. Lett. 68, 1160 (1992).
  * Wang  _et al._ (1998) S. C. Wang, U. Kurpick, and G. Ehrlich, Phys. Rev. Lett. 81 (1998).
  * Kellogg (1993) G. L. Kellogg, Appl. Surf. Sci. 67, 134 (1993).
  * Wen  _et al._ (1994) J. M. Wen, S.-L. Chang, J. W. Burnett, J. W. Evans, and P. A. Thiel, Phys. Rev. Lett. 73, 2591 (1994).
  * Pai  _et al._ (1997) W. W. Pai, A. K. Swan, Z. Zhang, and J. F. Wendelken, Phys. Rev. Lett. 79, 3210 (1997).
  * Giesen and Ibach (2003) M. Giesen and H. Ibach, Surf. Sci. 529, 135 (2003).
  * Giesen  _et al._ (1998) M. Giesen, G. S. Icking-Konert, and H. Ibach, Phys. Rev. Lett. 80, 552 (1998).
  * Fern  _et al._ (2000) J. Fern, L. Gmez, J. M. Gallego, J. Camarero, J. E. Prieto, V. Cros, A. L. V. de Parga, J. J. de Miguel, and R. Miranda, Surf. Sci. 459, 135 (2000).
  * der Vegt  _et al._ (1995) H. A. V. der Vegt, J. Alvarez, X. Torrelles, S. Ferrer, and E. Vlieg, Phys. Rev. B 52, 17443 (1995).
  * Busse  _et al._ (2003) C. Busse, C. Polop, M. Muller, K. Albe, U. Linke, and T. Michely, Phys. Rev. Lett. 91, 056103 (2003).
  * Muller  _et al._ (2005) M. Muller, K. Albe, C. Busse, A. Thoma, and T. Michelyi, Phys. Rev. B 71, 075407 (2005).
  * Liu and Adams (1992) C. L. Liu and J. B. Adams, Surf. Sci. 268, 73 (1992).
  * Tung and Graham (1980) R. T. Tung and W. R. Graham, Surf. Sci. 97, 73 (1980).
  * Fu and Tsong (2000) T. Y. Fu and T. T. Tsong, Surf. Sci. 454-456, 571 (2000).
  * Flahive and Graham (1980) P. G. Flahive and W. R. Graham, Surf. Sci. 91, 449 (1980).
  * Liu  _et al._ (1991) C. L. Liu, J. M. Cohen, J. B. Adams, and A. F. Voter, Surf. Sci. 253, 334 (1991).
  * Rice  _et al._ (1992) B. M. Rice, C. S. Murthy, and B. C. Garrett, Surf. Sci. 276, 226 (1992).
  * Stoltze (1994) P. Stoltze, J. Phys.: Condens. Matter 6, 9495 (1994).
  * Li and DePristo (1996) Y. Li and A. E. DePristo, Surf. Sci. 351, 189 (1996).
  * Mortensen  _et al._ (1996) J. J. Mortensen, B. Hammer, O. H. Nielsen, K. W. Jacobsen, and J. K. Norkov, “Elementary processes in excitations and reactions on solid surfaces,” (Springer-Verlag, Berlin, 1996) Chap. Density functional theory study of self-diffusion on the (111) surfaces of Ni, Pd, Pt, Cu, Ag and Au, pp. 173–182.
  * Chang  _et al._ (2000) C. M. Chang, C. M. Wei, and S. P. Chen, Phys. Rev. Lett. 85, 1044 (2000).
  * Haug and Jenkins (2000) K. Haug and T. Jenkins, J. Phys. Chem. B 104, 10017 (2000).
  * Kurpick (2001) U. Kurpick, Phys. Rev. B 64, 075418 (2001).
  * Bulou and Massobrio (2005) H. Bulou and C. Massobrio, Phys. Rev. B 205427 (2005).
  * Kim  _et al._ (2006) S. Y. Kim, I.-H. Lee, and S. Jun, Phys. Rev. B 74, 073412 (2006).
  * Repp  _et al._ (2003) J. Repp, G.Meyer, K. H. Rieder, and P. Hyldgaard, Phys. Rev. Lett. 91, 206102 (2003).
  * Trushin  _et al._ (2005) O. Trushin, A. Karim, A. Kara, and T. S. Rahman, Phys. Rev. B 72, 115401 (2005).
  * Karim  _et al._ (2006) A. Karim, A. N. Al-Rawi, A. Kara, T. S. Rahman, O. Trushin, and T. Ala-Nissila, Phys. Rev. B 73, 165411 (2006).
  * Nandipati  _et al._ (2009) G. Nandipati, Y. Shim, J. G. Amar, A. Karim, A. Kara, T. S. Rahman, and O. Trushin, J. Phys.: Condens. Matter 21, 084214 (2009).
  * Nandipati  _et al._ (2011) G. Nandipati, A. Kara, S. I. Shah, and T. S. Rahman, J. Phys.: Condens. Matter 23, 262001 (2011).
  * Shah  _et al._ (2012) S. I. Shah, G. Nandipati, A. Kara, and T. S. Rahman, J. Phys.: Condens. Matter 24, 354004 (2012).
  * Jónsson  _et al._ (1998) H. Jónsson, G. Mills, and K. W. Jacobsen,  _Classical and Quantum Dynamics in Condensed Phase Simulations_ , edited by B. J. Berne (World Scientific, Singapore, 1998).
  * Foiles  _et al._ (1986) S. M. Foiles, M. I. Baskes, and M. S. Daw, Phys. Rev. B 33, 7983 (1986).
  * Yildirim  _et al._ (2007) H. Yildirim, A. Kara, and T. S. Rahman, Phys. Rev. B 76, 165421 (2007).
  * Yildirim  _et al._ (2005) H. Yildirim, A. Kara, S. Durukanoglu, and T. S. Rahman, Surf. Sci. 600, 484 (2005).
  * Chirita  _et al._ (2000) V. Chirita, E. P. Munger, J. E. Greene, and J.-E. Sundgren, Thin Solid Films 370, 179 (2000).
  * Chirita  _et al._ (1999) V. Chirita, E. P. Munger, J. E. Greene, and J.-E. Sundgren, Surface Science 436, L641 (1999).
  * Einstein (1905) A. Einstein, Ann. Phys. (Leipzig) 17, 549 (English transl. Investigations on the Theory of Brownian Movement, Dover, New York, 1956) (1905).
  * Shah  _et al._ (shed) S. I. Shah, G. Nandipati, and T. S. Rahman, (to be published).


# Real-Time Joint Simulation of LiDAR Perception and Motion Planning for
Automated Driving

Zhanhong Huang, Xiao Zhang and Xinming Huang *This work was supported by the
US NSF Grant 2006738 and by The MathWorks. The authors are with the Department
of Electrical & Computer Engineering and Department of Computer Science,
Worcester Polytechnic Institute, Worcester, MA 01609, USA.
<EMAIL_ADDRESS>
###### Abstract

Real-time perception and motion planning are two crucial tasks for autonomous
driving. While there are many research works focused on improving the
performance of perception and motion planning individually, it is still not
clear how a perception error may adversely impact the motion planning results.
In this work, we propose a joint simulation framework with LiDAR-based
perception and motion planning for real-time automated driving. Taking the
sensor input from the CARLA simulator with additive noise, a LiDAR perception
system is designed to detect and track all surrounding vehicles and to provide
precise orientation and velocity information. Next, we introduce a new
collision bound representation that relaxes the communication cost between the
perception module and the motion planner. A novel collision checking algorithm
is implemented using line intersection checking that is more efficient for
long distance range in comparing to the traditional method of occupancy grid.
We evaluate the joint simulation framework in CARLA for urban driving
scenarios. Experiments show that our proposed automated driving system can
execute at 25 Hz, which meets the real-time requirement. The LiDAR perception
system has high accuracy within 20 meters when evaluated with the ground
truth. The motion planning results in consistent safe distance keeping when
tested in CARLA urban driving scenarios.

##  I INTRODUCTION

###  I-A Motivation

Perception and Motion Planning are two crucial tasks in an automated driving
system. The perception module in an autonomous vehicle (AV) receives the data
from the on-board sensors, process the information instantaneously, and send
the vehicle perception results to the motion planner. The task of autonomous
driving perception has been studied extensively using different approaches.
The traditional object based methods include object detection, object
classification and object tracking[1]. With advance of deep learning models,
frame based methods such as semantic segmentation[2], panoptic
segmentation[3], and 4D panoptic segmentation[4] have become popular. The
critical requirements of a perception system are real-time processing with
high accuracy and computational efficiency.

Upon receiving the information from the perception module, motion planner
makes behavior decisions in real-time, including path planning and velocity
planning, and then send commands to the executor, i.e., vehicle’s kinematic
controller and electronic control unit (ECU). Motion planner must be robust
and efficient. In order to meet the strict real-time requirement, it is
necessary to build a system in which the perception and motion planner modules
communicates and coordinates accordingly.

Figure 1: Illustration of our system pipeline. Bottom left: CARLA
simulator[5], Top Left: Lidar raw data shown in Rviz. Top Right: visualization
of the Lidar perception. Bottom Right: the visualization of the Motion
Planning

The existing works of perception were mainly focused on accuracy. Datasets
with metrics are defined as benchmarks for various perception tasks. For
object detection, a perception system must segment the objects and predict the
bounding box that aligns with the ground truth. Semantic segmentation and
panoptic segmentation utilize mIoU to evaluate the performance of point-wise
prediction results. Many perception research were aimed at improving the
benchmark metrics by a few percentages or even a fraction of one percent. The
question is - Does these progressive improvement on the perception benchmark
metrics affect the actual automated driving outcomes? In other words, does the
small increment in perception performance achieved by expensive computations
and energy cost really benefit the motion planning? This is the motivation of
our research work.

For trajectory planning, collision checking is often a time-consuming yet
indispensable step. A general approach in existing works requires discrete
representation of the environment in Frenet coordinate system[6, 7, 8, 9]. It
first uses the bounding circles to represent the collision constraint of all
vehicles and obstacles in the surrounding area on the map. Next, it transforms
the ego vehicle along the path on the map to check whether the obstacles would
appear in the ego vehicle’s bounding representation. The advantage of this
method is that the data structure of grid space representation is convenient
for path searching. Similar data structure can be directly applied to the
sensor input as a geometric representation of all objects on the road.
However, this method also has several disadvantages: Firstly, the cost of grid
representation is high which is directly proportional to the map resolution
and the sensing distance. Because of the high cost of mapping the environment,
areas outside of the driving lanes are usually neglected. However, objects not
on the ego lane, such as pedestrians and misbehaved vehicles at an
intersection, must be considered in the path planning process. Secondly, the
traditional method assumes that the geometric representation of the obstacles
all appears at once. But in practice the shape of an object may not reveal
completely in a LiDAR frame because of partial occlusion. Thirdly, the 2D grid
map is good for representation of positions, but not for other measurements
such as velocity. It would be time-consuming to convert all types of sensor
data into this grid mapping system.

Apart from the issues mentioned above, most existing works evaluates each
module in an automated driving system independently. For instance, the path
planning module always assumes the perception of the environment is perfect
and gets feedback instantaneously[10, 9]. The perception module is assumed
with a perfect localization, which is a lack of consideration since the two
modules are coupled in an AV system and affects each other interactively. Each
sub-module is optimized without being evaluating in a complete automated
driving framework. Hence, there is a possibility of leading to unexpected
behaviors when integrating these modules together during vehicle operations.

###  I-B Contribution

In contract to the existing works, the proposed framework validates the
integration of both perception and motion planning together in a real-time
autonomous driving simulation. The main contributions are summarized as
follows:

  * •

We build a real-time LiDAR perception system using point cloud clustering and
tracking-by-detection. An optimized non-iterative L-shaping method is designed
to estimate the vehicle position and orientation accurately. Furthermore, we
implement the vehicle velocity estimation via dynamic observations based on
the UKF-GNN tracking model.

  * •

We introduce a novel and efficient collision bound representation that reduces
the communication cost between the perception module and the motion planner.
The instance based collision checking can work independently without a grid
map representation.

  * •

We conduct real-time experiments in CARLA to validate the joint simulation
approach. The integrated system with LiDAR perception and motion planning
meets the urban driving performance requirement.

##  II Related Work

###  II-A Vehicle Detection and Tracking

Tracking-by-detection is a popular method for the task of multi-object
tracking (MOT). In this section we summarize the non-learning LiDAR-based
detection and tracking algorithms.

Detection of all surrounding objects is an important task in LiDAR point cloud
processing for data cleaning, candidate instancing, feature extraction, and
object classification. Voxel downsampling is one of the most popular ways to
pre-process the Lidar point cloud. The method preserves the spatial
distribution of the points while reducing the point size. Ground segmentation
is another essential step in pre-processing. RANSAC plan fitting [11] and
slope-based segmentation proposed by [12] are two of the representative
techniques. They describe the ground as a single plane or continuously
extended surface in vertical and horizontal directions. Subsequently, a
traditional method including instancing, feature extraction, and object
classification is often implemented step-by-step. Point cloud clustering is
proved to be an efficient way to instance the object candidates in the point
set. Distance-based, density-based, and range image-based methods like
Euclidean clustering, DBSCAN [13] and depth clustering[14] can be selectively
applied for Lidar data[15]. The candidates are then processed by the
classifier to tag semantic labels of the objects. Through the recent works[16]
[17], support vector machine (SVM) [18], random forest (RF)[19] and simple
rule-based method can efficiently classify the point cloud clusters.

Tracking is a step that takes the detected objects, tracking and identifying
them as unique samples using the historical information sequence. This is done
through a framework that combines object state prediction (filter) and data
association to deal with object dynamics. The extended Kalman filter (EKF) and
Unscented Kalman Filter (UKF) are commonly used nonlinear approaches for
filtering, with the UKF outperforming the EKF in predicting heavy nonlinear
cases.[20]. Other filters, such as interactive multiple model (IMM) filter and
particle filter can also be used as alternative solutions. [21]. The Lidar-
based vehicle tracking studies typically implemented the global nearest
neighbor (GNN) and joint probability distribution association filter (JPDAF)
methods for data association. The JPDAF outperforms the GNN in complex cases
where multiple measurements are very close [17].

###  II-B Collision Checking in Trajectory Planning

The existing collision checking methods for automated driving can be
categorized into two main categories, i.e. footprint based and potential field
based.

Footprint based methods use a few circles along the ego vehicle axis to bound
the vehicle to create the collision bounding, where the number of circles and
and their radii vary in different representations. Zigler and Stiller[22]
introduced a fast collision checking method, where they decomposed the vehicle
shape into several disks and then converted the disks into a rectangle shape
to fit in the grid map. Zhang et al.[6] proposed a hybrid trajectory planning,
where they evaluated the clearance of vehicle footprint in the grid map, which
was computed based on the distance between the way point and the obstacle in
each specific direction. This clearance allows the path optimizer to find the
orientation of the path since the cost increases when the clearance gets
shorter.

Potential field based methods calculated an ellipse around each object, where
the ellipse surrounds the reachable pose in several future time steps provided
with the object’s current velocity and acceleration. Philipp and Goehring[23]
proposed a collision octagon check boundary that facilitates the calculation
of collision state probability and collision event probability. The octagon
representation simplified the bounding rectangle of two possible colliding
vehicles to a point and an octagon, where the octagon’s vertex was calculated
from the Gaussian state distribution of the obstacle.

##  III Joint Simulation Framework

Figure 2: Our joint simulation system is based on the CARLA framework and ROS
bridge. In the proposed system, the Lidar perception and the motion planner
are realized as ROS nodes and work asynchronously. The dashed lines show the
asynchronous signals among ROS nodes and the solid lines indicate synchronous
signals with each node.

Fig. 1 shows our system pipeline and Fig. 2 outlines the proposed system
framework with all modules and the connections among them. The perception
module takes the converted CARLA sensor messages from ROS Bridge as input[5].
It generates the environment abstractions that allow the decision module in
the local planner to decide whether it needs to request global waypoints from
the global planner. The module’s decision is usually based on the road
policies, such as road speed limit changes or ego vehicles stopped by traffic
lights. The path generation process starts once the global waypoints are
retrieved from the global planner. The path generation module first requests
the surrounding obstacles from the Lidar perception module and performs the
collision checking. Subsequently, it outputs the collision-free path to the
velocity planner for speed estimation and the decision module for loop closure
checking. Next, the target pose and speed is sent to the controller to
generate the control output. Finally, the control output is sent to the ROS
bridge for vehicle movement and display update in CARLA.

##  IV Technical Approach

In this section, we first present the proposed Lidar perception module and its
algorithms. Then, we introduce the motion planner module, which includes the
collision checking algorithm and the velocity model.

###  IV-A System Assumptions

For the perception system, we assume the ego car’s speed and angular velocity
can be fetched by car-deployed sensors such as the wheel tachometer and IMU.
Traffic signal status and stop signs can be retrieved from the the image-based
sensors or V2X communications. For the motion planner, we assume our global
waypoints have already been calculated by a remote server, which is not
required to be handled by the ego vehicle. We further assume that the ego
vehicle gets its global position through GPS, and its position estimation is
more accurate than the waypoints resolution.

###  IV-B Lidar Perception

####  IV-B1 Pre-processing

As the front-end of a real-time and high-performance Lidar perception system,
pre-processing is applied to remove the noise and provide computationally
efficient point clouds for downstream processing.

Down-sampling. The time complexity of point cloud algorithms largely depends
on the quantity of points. A standard Lidar frame with 64 channels typically
consists of approximately 20,000 to 30,000 points, which can be
computationally intensive if all points are processed.

Depending on the applications, not all points are equally important. To focus
on the potentially moving objects near the ego vehicle, the pre-processing
step segments out an origin-centered region of interest (RoI) from a bird-eye
view (BEV). Next, a voxel down-sampling method is applied to reduce the
negative effects of outliers on the subsequent segmentation and classification
processes.

Ground Segmentation. As the ground points are naturally connected to nearly
all objects in the Lidar scan, distance or density-based point cloud
instancing methods do not perform well before ground removal. To address the
problem and further reduce the number of points, the RANSAC plane fitting
technique is utilized to separate the potential vehicle candidates and the
ground plane. Although this iterative RANSAC algorithm is not particularly
time efficient, ground plane removal reduces the number of points, which
speeds up the subsequent processing considerably.

####  IV-B2 Vehicle Detection

We propose an efficient object detection algorithm based on point cloud
clustering that does not require deep learning. It effectively extracts the
vehicle candidates in each frame.

Vehicle Instancing. To identify the vehicle candidates from non-ground Lidar
points, we implement the DBSCAN method to cluster points based on their
spatial density distribution. By configuring the density parameter $\epsilon$
and minimum points within the density region, DBSCAN extracts the clusters
with a specific density in a certain range and other sparser distributed
points are labeled as noise.

L-shape Fitting. From our study, the essential requirements of a perception
system in local planning for autonomous driving involve determining the
position and orientation of the detected objects. To extract the information
from each cluster of points, we adopt the L-shape fitting approach. The method
fits a rectangle along the L edge from BEV. While L-shape fitting via variance
is suitable for 2D point clouds and offers optimal vehicle box alignment, it
cannot be directly applied to 3D points due to the large number of points
which would require much longer computing time. Instead, we project each
cluster of points to BEV and then extract the convex shape of the points. This
technique optimizes the run time while preserving the L-shape fitting
performance.

Rule-based Classifier. Rule-based classifier is an efficient method for
filtering out vehicle candidates based on the geometric and statistical
criteria. For vehicle classification, it is highly efficient and provides
reliable results when rule criteria are specified accordingly. In this work,
we set the thresholds for various parameters such as the number of points,
width, length, area, and the width-to-length ratio of the bounding box from
BEV for vehicle detection.

###  IV-C Multi-Vehicle Tracking

UKF-GNN Framework. Multi-Vehicle tracking is similar to the classic task of
multi-object tracking. The traditional solution comprises two key components:
object behavior predictor and data association framework. We implemented the
UKF as the predictor while using GNN as the data associator. When the new
candidates sequences are fed, the data associator calculates the correlation
scores among the objects in the detected list and predicts objects in the
tracking list. Based on the score matrix, the tracking list is updated
according to the mapping and list management policy. During the update step,
the Kalman filter of each object takes the matched detection object
information as measurement to update its probability of the states model and
subsequently predicts the location and speed of the object for the next state.

###  IV-D Dynamic Velocity Estimation.

Estimating the velocity of objects from a moving vehicle is considerably more
complex than the estimation from stationary observations. This is due to the
fact that the current observation coordinates undergo both rotation and
translation in the 3D space to its previous observation. Therefore, the
initial step in states updating is to adjust the previous observation states
using the ego car’s speed and angular measurement described by the following
equations:

$R(\theta,\Delta x,\Delta y)=\begin{bmatrix}\cos{\theta}&-\sin{\theta}&\Delta
x\\\ \sin{\theta}&\cos{\theta}&\Delta y\\\ 0&0&1\end{bmatrix}$ (1)

$d_{ego}=v_{ego}\Delta t$ (2)

$\begin{split}\Delta x&=-\sin{\theta}d_{ego}\\\ \Delta
y&=-\cos{\theta}d_{ego}\end{split}$ (3)

$\begin{bmatrix}x_{t}\\\ y_{t}\\\ 1\end{bmatrix}=R(\theta_{ego},\Delta
x,\Delta y)\begin{bmatrix}x_{t-1}\\\ y_{t-1}\\\ 1\end{bmatrix}$ (4)

where $\theta\subseteq[\pi,-\pi)$ is the orientation of ego car in the map
coordinate, $\Delta t$ is time difference between each frame. For each lidar
frame, we have recorded the position of the detected vehicles. Once the object
positions in the previous frame are transformed to the current one, we can
estimate the velocities and orientations of the tracked vehicles accurately.

###  IV-E Motion Planner

####  IV-E1 Waypoint conversion and optimization

Before applying the perception results for collision checking, we first need
to convert our waypoints retrieved from the global path planner into the ego
car frame. We consider the path starting from the frame origin, where the
origin is set on the ego car’s ground center. The frame coordinates are mapped
with the right-hand rule, and the _x_ positive direction is set to the
orientation of the ego car. Given the coordinate frame of the ego car, we can
map these waypoints into the ego frame with the transformation matrix:

$\begin{bmatrix}wp_{eg,x}\\\ wp_{eg,y}\\\
1\end{bmatrix}=R(\theta_{ego},-p_{x},-p_{y})\begin{bmatrix}wp_{mp,x}\\\
wp_{mp,y}\\\ 1\end{bmatrix}$ (5)

where $\theta_{ego}\subseteq[\pi,-\pi)$ is the orientation of the ego car in
the map frame, and $p$ is the position of the ego car on the map. So that we
can get a series of waypoints located in the ego frame:

$wp=[wp_{0},\dots,wp_{m}]$ (6)

where $m$ is the minimum distance for perception, calculated as the safe
braking distance based on the ego car’s speed:

$d_{brake}={{v_{ego}}\over{2\mu g}}$ (7)

where $\mu$ is the friction coefficient, which we set as a constant equal to
0.35 in our simulation. _g_ is the gravitational acceleration, which we also
set as a constant equal to 9.8 m/$\text{s}^{2}$ in our simulation. Given the
safe braking distance, we can calculate the number of converted waypoint in
each cycle by:

$m={{d}\over{d_{wp}}}\times f_{safe},\;m\subseteq\mathbf{N}$ (8)

where the $d_{wp}$ is the preset distance between two waypoints, we applied
the safety factor $f_{safe}$ to get an extra safety buffer.

Once we get the converted waypoints, we need to use them to calculate a smooth
trajectory. In our design, we implement B-Spline to generate the 2-degree
Bézier curves, in which the waypoints’ _x_ and _y_ positions are being
considered separately, and the reference axis is the travel euclidean
distance. Given the reference axis, we can calculate the B-Spline for our
trajectory, referring to the x and y axes. Then we can re-sampled the x and y
trajectory to get our re-sampled waypoints given a fixed travel distance
increment $t_{s}$, which is set to 0.5 in our simulation:

$RWpList=[(trajectory_{x}(t_{s}),trajectory_{y}(t_{s})]$ (9)

####  IV-E2 Collision checking

Next, our motion planner performs collision checking along the trajectory. In
our design, we make an intersection checking between our trajectory and the
car bounding boxes. Compared to traditional bubble collision checking, our
algorithm is more computationally efficient for long range planning. Moreover,
the simplified representation of obstacles can be directly obtained from the
lidar perception, which does not need to convert into the BEV. Finally, line
intersection checking is more utilized in the detected object motion
representation, where the trajectory of the dynamic obstacle and the obstacle
itself can be evaluated together.

Based on our assumption that an observed vehicle always drives toward its
heading, we can represent the trajectory of a vehicle by extending the front
bounding points, which can be calculated using its speed and orientation:

$\begin{split}E_{x}&=\cos{\theta_{c}}speed_{c}t_{est}\\\
E_{y}&=\sin{\theta_{c}}speed_{c}t_{est}\end{split}$ (10)

where $theta_{c}$ is the vehicle orientation in the ego frame, $speed_{c}$ is
the vehicle current speed, and $t_{est}$ is the estimated time for the vehicle
trajectory. An example of the vehicle bounding extension model is shown in
Fig. 3.

Figure 3: An example of the vehicle velocity extended bounding model. The red
lines are the lidar detected bounding and the orange lines are the velocity
extended bounding. The green lines are the ground truth bounding extended by
the estimated velocity.

By adding the extended distance to the front-end bounding points, we can get
an extended bounding box that includes the velocity information in the current
time frame. Then we can generate Bounding Lines around the vehicle through the
bounding coordinates:

$\begin{split}VBound&=[L_{1},L_{2},L_{3},L_{4}]\\\ &=[L_{c1\to c2},L_{c1\to
c3},L_{c2\to c4},L_{c3\to c4}]\end{split}$ (11)

Figure 4: Frome left to right: (a) Original global waypoints. (b) Resampled
Waypoints after Bspline, trajectory line with side bound assisted check lines
(c) Retrieved the interesect point with lowest index. (d) Collision-free path.

Next, we sort the $VBound$ along with the information of each vehicle, which
is also shown in lines one to six of Algorithm 1. Then we perform the
intersection check along the trajectory, which is introduced in line seven of
Algorithm 1 and detailed described in Algorithm 2 also shown in Fig. 4.
Eventually, after the collision checking, we get the distance, speed, and
orientation of the vehicle lying on our trajectory. A new collision-free
trajectory is generated.

Algorithm 1 Generate Collision Free Path with vehicle info

1:$VList,RWpList$

2:$distance,speed,\theta,newPath$

3:

4:Init $BoundingList$

5:Init $InfoList$

6:for each $Vehicle$ in $VList$ do
$\triangleright$ Retrieve vehicle information and bounding box line

7: Add List[$bounds$] to $BoundingList$

8: Add List[$speed,yaw,position$] to $InfoList$

9:end for

10:

11:$VIndex,distance,newPath\newline
\leftarrow\textsc{IntersectCheck}(RWpList,BoundingList)$
$\triangleright$ Get closest vehicle index that occupy on our path, its
distance to us and our collision free path

12:

13:$\theta\leftarrow InfoList[VIndex][1]$

14:$speed\leftarrow InfoList[VIndex][2]$
$\triangleright$ Retrieve vehicle measurement from info list

15:

16:return $distance,speed,\theta,newPath$

Algorithm 2 IntersectCheck

1:$RWpList,BoundingList$

2:$VIndex,CFWpList$

3:

4:Init $CFWpList\leftarrow RWpList$

5:for each $VBoundSet$ in $BoundingList$ do

6: for each $Bound$ in $VBoundSet$ do

7: Init $CheckFlag\leftarrow$ True

8: while $CheckFlag\;$and $CFWpList$ do

9: $Lines$ generate from $CFWpList$

10: if $Bound$ Intersect $Lines$ then

11: $CFWpList$.PopRight()

12: $VIndex\leftarrow$Index($VBoundSet$)

13: $distance\leftarrow$D($Bound,RWPList[0]$)

14: else

15: $CheckFlag\leftarrow$ False

16: end if

17: end while

18: end for

19:end for

20:return $VIndex$, $distance$, $CFWpList$

####  IV-E3 Speed Planning

After obtaining the collision-free path and the occupied vehicle information,
we need to determine the ego car speed and target pose as the input for the
PID controller. In the scenario of lane keeping, given the target pose, speed
is calculated for the following three cases:

$v_{pre}=\begin{cases}v_{n}&\text{no obstacle}\\\ v_{obs}&\text{obstacle in
front}\\\ v_{platoon}&\text{following vehicle}\end{cases}\\\ $ (12)

When there is no obstacle and vehicle in front of the ego car, we calculate
the speed by subtracting the estimated reachable distance and the target pose
distance multiplied by control timestep $\Delta t$:

$v_{n}=v_{current}+(d_{pose}-d_{reach})*\Delta t\\\ $ (13)

When the Intersection check algorithm reports an obstacle in front of the ego
car or a vehicle not in the same orientation as the ego car, we cap the speed
by estimating the difference between the distance from the obstacle and safe
distance $d_{safe}$:

$v_{obs}=v_{appr}+{d_{obs}-d_{safe}\over{d_{safe}}}v_{appr}$ (14)

where $v_{appr}$ is the preset speed for slowing down the ego car. And the
safe distance $d_{safe}$ is calculated by $d_{brake}$ in equation 7, and an
additional buffer distance $d_{buffer}$ to keep the minimum distance:

$d_{safe}=d_{brake}+d_{buffer}$ (15)

When the intersection algorithm finds that there is a leading vehicle in front
of the ego car, our target speed is calculated by

$v_{platoon}=v_{lead}+w(d_{lead}-d_{safe})\Delta t$ (16)

where $v_{lead}$ is the leading vehicle speed, $d_{lead}$ is the distance
between the leading vehicle and the ego car, and $w$ is a proportional weight
to apply to the distance difference. Once the preferred speed is calculated,
we bound the speed and send it to the controller.

$v_{exc}=\begin{cases}0&v_{pre}<0\\\ v_{pre}&0\leq v_{pre}<v_{reach}\\\
v_{reach}&v_{pre}\geq v_{reach}\end{cases}$ (17)

where the maximum reachable speed in the next control timestep $v_{reach}$ is
calculated with

$v_{reach}=v_{current}+a_{max}*\Delta
t,\;v_{reach}\subseteq[v_{init},v_{max}]$ (18)

where $v_{max}$ is the speed limit on the road, and the $v_{init}$ is the
initial launch speed that allows the the car ramp up the speed faster from
zero.

##  V Empirical Results

###  V-A Implementation Setup

We first initialized the CARLA simulator and the CARLA ROS Bridge, which
handles the connection between our system and CARLA[5, 24]. We choose Town10
and Town01 as the simulation maps in CARLA with 50 vehicles but no
pedestrians. We set a fixed 50 ms as our simulation time-step, ticked by the
ROS bridge. For the vehicle speed limit, we set it to 8.33 m/s for Town10, and
adaptive max speed based on the speed limit posted on the signs as in Town01.
For the $v_{init}$, the initial launch speed, we set it to 6 m/s. The global
waypoints were generated by CARLA default waypoint publisher with 2 m
resolution. For the trajectory optimization, we re-sampled the trajectory with
0.5 m as the fixed travel distance step $t_{s}$. For the speed planning, we
set the ego vehicle parameter to $\mu=0.35$, $a_{max}$ = 2.5 m/$\text{s}^{2}$.
We use the CARLA default PID controller which is already tuned with the ego
vehicle as a Tesla Model 3. The proposed method is implemented in Python,
where the lidar perception module implements Open3d library without GPU
acceleration. Each module ran in a ROS1 node and communicated with each other
with ROS service protocol. All the experiments were conducted on a desktop
computer with an AMD 5950X CPU (4.5GHz).

###  V-B Motion Planner Verification

Figure 5: Distance to the front vehicle during automated platooning.

To validate our motion planner performance, we let the ego vehicle run in both
maps with random global waypoints for 40 minutes using ground truth dynamic
objects information provided by CARLA. We recorded the distance between ego
vehicle and the detected front vehicle in platooning scenarios, which is shown
in the right-hand side of Fig. 5. We can see that the main distribution of the
recorded following distance in Town01 is longer than that of Town10. This is
because the speed limit in Town01 is adaptive, which can reach a maximum speed
of 15 m/s. For this reason, the closest distance in Town01 is shorter than in
Town10.

We also recorded the execution time for path generation and velocity
calculation for the reference of system evaluation in the later experiments,
shown in the first two columns of Table.I. We can see the whole system running
above 20 Hz. The most time-consuming step is to transform the global waypoints
to the ego frame and then back to the global frame after processing. The cases
of long processing time occurred when there were too many resampled waypoints
to be converted. Because the waypoint number is directly related to the
vehicle speed, we can also see the average processing time of Town01 is less
than that of Town10.

TABLE I: Time consumption for each step |  | Town01 GT |  | Town10 GT |  | Town01 Lidar |  | Town10 Lidar
---|---|---|---|---|---|---|---|---
Time |  | Avg. (s) | Max. (s) |  | Avg. (s) | Max. (s) |  | Avg. (s) | Max. (s) |  | Avg. (s) | Max. (s)
Lidar Perception |  | / | / |  | / | / |  | 0.041 | 0.101 |  | 0.034 | 0.090
Path Generation |  | 0.009 | 0.045 |  | 0.012 | 0.053 |  | 0.009 | 0.027 |  | 0.011 | 0.024
Speed Planning |  | 0.001 | 0.006 |  | 0.002 | 0.03 |  | 0.001 | 0.008 |  | 0.001 | 0.003
Local Planner Total |  | 0.044 | 0.120 |  | 0.028 | 0.076 |  | 0.050 | 0.143 |  | 0.031 | 0.092

###  V-C Lidar Perception Evaluation

We evaluate our perception by three different kinds of metrics. For each
experiment, the data is collected by more than 5000 consecutive lidar frames
in Town01. In Tab. II, $Recall$ indicates the truth positive detection of the
surrounding vehicles by the LiDAR perception system. As the effective range
decreases, the accuracy of the perception system becomes higher, thus the
safety of the autonomous driving system is better. The 20-meter perception
range setup is adapted according to the maximum speed and minimum braking
distance. The $IoU$ represents the ratio of the intersection area between the
detected and ground truth bounding box over the union area. The first two rows
of overall object detection are relatively lower than the last two rows of
moving object detection since we implement different strategies to align the
bounding box to dynamic and static (speed estimated around zero) objects. From
$E_{\theta}$, our optimized L-shape detection algorithm has a very high
accuracy for estimating the object’s orientation, with only a 2-3 degrees of
estimation error. Speed estimation error $E_{v}$ is also under 0.6 $m/s$. When
the perception system runs at 20 $Hz$, each time step is 0.05 second and the
speed error is possibly caused by the LiDAR measurement noise.

TABLE II: Evaluation for Lidar Perception

Position Orientation (deg) Speed (m/s) Range $Recall$ $mIoU$ $E_{\theta}$
(Avg./Std.) $E_{v}$ (Avg./ Std.) $R$ = 20m 0.83 0.62 2.03 / 0.11 0.49 / 0.85
$R$ = 15m 0.97 0.64 2.12 / 2.18 0.41 / 0.61 Dyn. $R$ = 20m 0.94 0.71 3.08 /
6.05 0.59 / 2.23 Dyn. $R$ = 15m 0.99 0.73 1.57 / 3.77 0.55 / 0.82

###  V-D System Evaluation

To evaluate the system performance by combing the lidar perception module with
the motion planner, we ran the same simulations with the ground truth
perception and recorded the platooning following distance and total execution
time. The following distances are shown in the left-hand side of Fig.5, and
the time consumption of each step is also shown in Table.I. The overall time
consumption is only increased by around 6 ms compared to the time consumption
with ground truth information thanks to the asynchronous design for each node
running in the perception framework.

The run-time table shows that lidar perception processing has an apparent
slowdown compared to Town10 to Town01. The result is caused by the building
and environmental features. Town01 is a country town which includes many
little infrastructures like fences. It leads to over segments of vehicle
candidate instancing as well as a rise in burden to the downstream processing
like box alignment and tracking data association. But owing to the
conservative policy design, the perception system can still achieve a zero-
collision record in the 40 minutes long simulations.

We also compared the platooning performance between different simulations
within the same map. On the top of Fig. 6 shows one example of the system with
ground truth information. And the bottom of Fig. 6 shows one case with the
Lidar perception. In both figures, the red, blue, and green lines represent
the ego vehicle speed, the following vehicle speed, and the following
distance. Compared two scenarios, we can see both systems performed a relative
distance keeping as the speed increased. But the speed estimation from the
Lidar perception module had more noise than the ground truth information,
which caused variations in the distance keeping.

Figure 6: Top: Following distance data example by system using ground truth
information. Bottom: Following distance data example by system with Lidar
Perception.

##  VI Conclusion

In this paper, we present an autonomous vehicle joint simulation system with
both lidar perception and planning module deployed. Firstly, a traditional
tracking-by-detection framework is employed to handle the lidar-based object
detection. We optimize the L-shape detection algorithm and simplified the
perception pipeline. Next, we introduce a new collision checking algorithm
that allows relaxation of obstacle representations while also improves the
long range collision check efficiency. We evaluate our system performance by
simulating urban driving scenarios in CARLA simulator. Similar to the
simulations using ground truth vehicle data, our system with LiDAR perception
can run at an average of 20 Hz that meets the real-time requirement. Finally,
within the range of 20 meters, the perception system gives accurate vehicle
positions with a recall above 0.83 and mIOU above 0.62. Based on the dynamic
velocity estimation, the average error of orientation is less than 3.08 degree
and the average speed estimation error is less than 0.59 m/s.

## References

  * [1] A. Geiger, P. Lenz, and R. Urtasun, “Are we ready for autonomous driving? the kitti vision benchmark suite,” in _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2012.
  * [2] J. Behley, M. Garbade, A. Milioto, J. Quenzel, S. Behnke, C. Stachniss, and J. Gall, “SemanticKITTI: A Dataset for Semantic Scene Understanding of LiDAR Sequences,” in _Proc. of the IEEE/CVF International Conf. on Computer Vision (ICCV)_ , 2019.
  * [3] J. Behley, A. Milioto, and C. Stachniss, “A benchmark for lidar-based panoptic segmentation based on kitti,” _arXiv preprint arXiv:2003.02371_ , 2020.
  * [4] M. Aygün, A. Osep, M. Weber, M. Maximov, C. Stachniss, J. Behley, and L. Leal-Taixe, “4D Panoptic Segmentation,” in _In Proc. of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2021.
  * [5] A. Dosovitskiy, G. Ros, F. Codevilla, A. Lopez, and V. Koltun, “CARLA: An Open Urban Driving Simulator,” Nov. 2017, arXiv:1711.03938 [cs]. [Online]. Available: http://arxiv.org/abs/1711.03938
  * [6] Y. Zhang, H. Chen, S. L. Waslander, J. Gong, G. Xiong, T. Yang, and K. Liu, “Hybrid Trajectory Planning for Autonomous Driving in Highly Constrained Environments,” _IEEE Access_ , vol. 6, pp. 32 800–32 819, 2018, conference Name: IEEE Access.
  * [7] L. Zheng, P. Zeng, W. Yang, Y. Li, and Z. Zhan, “Bézier curve-based trajectory planning for autonomous vehicles with collision avoidance,” _IET Intelligent Transport Systems_ , vol. 14, no. 13, pp. 1882–1891, 2020, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1049/iet-its.2020.0355. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1049/iet-its.2020.0355
  * [8] W. Xu, Q. Wang, and J. M. Dolan, “Autonomous Vehicle Motion Planning via Recurrent Spline Optimization,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , May 2021, pp. 7730–7736, iSSN: 2577-087X.
  * [9] C. Jie, C. Yingbing, Z. Qingwen, G. Lu, and L. Ming, “Real-Time Trajectory Planning for Autonomous Driving with Gaussian Process and Incremental Refinement,” May 2022, arXiv:2205.11853 [cs]. [Online]. Available: http://arxiv.org/abs/2205.11853
  * [10] R. Xu, H. Xiang, X. Han, X. Xia, Z. Meng, C.-J. Chen, C. Correa-Jullian, and J. Ma, “The OpenCDA Open-source Ecosystem for Cooperative Driving Automation Research,” _IEEE Transactions on Intelligent Vehicles_ , pp. 1–13, 2023, conference Name: IEEE Transactions on Intelligent Vehicles.
  * [11] B. Wang, J. Lan, and J. Gao, “Lidar filtering in 3d object detection based on improved ransac,” _Remote Sensing_ , vol. 14, no. 9, 2022. [Online]. Available: https://www.mdpi.com/2072-4292/14/9/2110
  * [12] P. Narksri, E. Takeuchi, Y. Ninomiya, Y. Morales, N. Akai, and N. Kawaguchi, “A slope-robust cascaded ground segmentation in 3d point cloud for autonomous vehicles,” in _2018 21st International Conference on Intelligent Transportation Systems (ITSC)_ , 2018, pp. 497–504.
  * [13] E. Schubert, J. Sander, M. Ester, H. P. Kriegel, and X. Xu, “Dbscan revisited, revisited: why and how you should (still) use dbscan,” _ACM Transactions on Database Systems (TODS)_ , vol. 42, no. 3, pp. 1–21, 2017.
  * [14] I. Bogoslavskyi and C. Stachniss, “Fast range image-based segmentation of sparse 3d laser scans for online operation,” in _2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2016, pp. 163–169.
  * [15] Y. Zhao, X. Zhang, and X. Huang, “A technical survey and evaluation of traditional point cloud clustering methods for lidar panoptic segmentation,” _arXiv preprint arXiv:2108.09522_ , 2021.
  * [16] J. Zhang, W. Xiao, B. Coifman, and J. P. Mills, “Vehicle tracking and speed estimation from roadside lidar,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 13, p. 5597–5608, 2020, conference Name: IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing.
  * [17] M. Sualeh and G.-W. Kim, “Dynamic multi-lidar based multiple object detection and tracking,” _Sensors_ , vol. 19, no. 6, p. 1474, Jan 2019, number: 6 Publisher: Multidisciplinary Digital Publishing Institute. [Online]. Available: https://www.mdpi.com/1424-8220/19/6/1474
  * [18] Z. Yan, T. Duckett, and N. Bellotto, “Online learning for human classification in 3d lidar-based tracking,” in _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , Sep 2017, p. 864–871.
  * [19] B. Lv, H. Xu, J. Wu, Y. Tian, Y. Zhang, Y. Zheng, C. Yuan, and S. Tian, “Lidar-enhanced connected infrastructures sensing and broadcasting high-resolution traffic information serving smart cities,” _IEEE Access_ , vol. 7, pp. 79 895–79 907, 2019.
  * [20] A. Arya Senna Abdul Rachman, “3d-lidar multi object tracking for autonomous driving: multi-target detection and tracking under urban road uncertainties,” 2017.
  * [21] T.-D. Vu and O. Aycard, “Laser-based detection and tracking moving objects using data-driven markov chain monte carlo,” in _2009 IEEE International Conference on Robotics and Automation_ , 2009, pp. 3800–3806.
  * [22] J. Ziegler and C. Stiller, “Fast collision checking for intelligent vehicle motion planning,” in _2010 IEEE Intelligent Vehicles Symposium_ , Jun. 2010, pp. 518–522, iSSN: 1931-0587.
  * [23] A. Philipp and D. Goehring, “Analytic Collision Risk Calculation for Autonomous Vehicle Navigation,” in _2019 International Conference on Robotics and Automation (ICRA)_ , May 2019, pp. 1744–1750, iSSN: 2577-087X.
  * [24] Stanford Artificial Intelligence Laboratory et al., “Robotic operating system.” [Online]. Available: https://www.ros.org