# REFERENCES

Clusters of interstitial carbon atoms

near the graphite surface as a possible origin

of dome-like features observed by STM

V.F. Elesin and L.A. Openov

Moscow State Engineering Physics Institute (Technical University),

Moscow, 115409, Russia

Abstract

> Formation of clusters of interstitial carbon atoms between the surface and second atomic layers of graphite is demonstrated by means of molecular dynamics simulations. It is shown that interstitial clusters result in the dome-like surface features that may be associated with some of the hillocks observed by STM on the irradiated graphite surface.

1. Introduction

Ion irradiation of graphite surface results in formation of defects that are seen by scanning tunneling microscopy (STM) . A single-ion impact creates a hillock whose nature still remains controversial since STM probes the electronic states of the surface rather than the actual surface topography. So, it is difficult to deduce the atomic structure in the vicinity of defects based on STM data alone.

A number of suggestions have been made concerning the interpretation of hillocks seen by STM on graphite surfaces irradiated by low-energy ($`100`$ eV) ions, including, e.g., the incident ions trapped between the surface and second layer of atoms in graphite , the vacancies on the graphite surface , self-atom implantation , etc. It seems, however, that in many cases the formation of hillocks involves carbon atoms only.

Recent molecular dynamics simulations  indicate that carbon interlayer interstitials below the surface layer account for $``$ 50 $`\%`$ of all defects produced by incident and secondary recoils with the energy $`E=300`$ eV $`÷`$ 3 keV, while $`100`$ $`\%`$ of defects produced by secondary recoils with the energy $`E=100`$ eV $`÷`$ 300 eV are single carbon interstitials. It was noticed in  that the single interstitials may migrate to form clusters between the layers, which may be the source of some of the hillocks. The authors of  arrived at a conclusion that, at most, half of the surface defects may be formed by interstitial carbon clusters.

As to high-energy ($`100`$ keV) ion implantation, the hillock shape of the defects was attributed by several authors to the formation of single carbon interstitials and interstitial carbon clusters between the graphite layers, each hillock being originated from a single ion impact . Two types of hillocks were found in , one of which showing the undisturbed (i.e. definitely without vacancies) lateral atomic arrangement on the (0001) graphite surface, while the other one exhibited a distorted atomic structure (but also presumably without vacancies). So, it is likely that the hillocks formed under both low- and high-energy ion irradiation of graphite are mainly due to interlayer carbon interstitials and/or their clusters. If this explanation is true, a unique opportunity appears to study the processes of cluster formation and a related phenomenon of radiation-induced swelling of graphite by means of STM. Besides, there is a need for a theory of formation, structure, and arrangement of such clusters. The theory should give an answer to several key questions: 1) what is the reason for attraction of interlayer interstitials; 2) what are the physical conditions of cluster formation; 3) why the recombination of interstitials with vacancies does not inhibit the formation of clusters.

The purpose of this paper is to study numerically the formation of clusters of interstitial carbon atoms in the interlayer region between the surface and second graphite layers. Recently we have shown by means of molecular dynamics simulations that interlayer carbon interstitials do form clusters in the bulk of graphite . The physical reason for this effect  is the deformation interaction of interstitials mediated by the nearby graphite sheets. An interstitial stretches the lattice in its vicinity, and the stretched domain appears to be attractive for other interstitials located within the same interlayer region, see Fig.1. In other words, the formation of interstitial clusters minimizes the deformation energy relatively to the case of spatially separated interstitials due to minimization of the overall buckling of adjacent graphite sheets. It is important that mobility of interstitial clusters is much lower than that of single interstitials, and hence the recombination of interstitials with vacancies is strongly suppressed .

On general grounds, one can expect the same mechanism of clusters formation to be realized near the graphite surface. Indeed, the formation of clusters of noble-gas atoms inserted between the topmost and second atomic layers of graphite has been numerically demonstrated by Marton et al. . We are not aware of similar studies on the clusters of carbon atoms.

The paper is organized as follows. First we present the general theory of diffusion instability of the uniform distribution of defects in solids. Next we numerically simulate the formation of interstitial carbon clusters underneath the surface layer of graphite. Our calculations show that such clusters can indeed be formed. Their specific shapes are dictated by the initial distribution of interstitials in the interlayer region and by the strong covalent interaction between the interstitials. The interstitial clusters give rise to dome-like features on the graphite surface that are similar to those observed experimentally by STM.

2. Diffusion instability of a system of defects in solids

In this Section we sketch the main features of the theory of diffusion instability that can lead to formation of clusters of defects in disordered solids, see  for more details. We shall demonstrate that the uniform distribution of defects is unstable with respect to separation in defect-rich and defect-poor regions if the concentration of defects is sufficiently high and exceeds some critical value that depends on material parameters and temperature. The mechanism of instability is rather general: the deformation of the crystal lattice by defects results in anomalous flow of defects along the gradient of defect concentration $`n(𝐫,t)`$, i.e. opposite in direction to the diffusion flow of defects.

In the simplest case of defects of the same type, the kinetic equation for $`n(𝐫,t)`$ is

$$\frac{n}{t}=QR(n)div𝐣,$$
(1)   
where $`Q`$ is the source of defects, the term $`R(n)`$ accounts for recombination processes, and $`𝐣`$ is the defect current density. In the case of isotropic continuum media, the expression for $`𝐣`$ reads

$$𝐣=D\frac{n}{𝐫}+\frac{D\mathrm{\Omega }n}{3T}\frac{\sigma _{ii}}{𝐫},$$
(2)   
where the first term corresponds to the diffusion of defects, while the second one describes the motion of defects in the field of elastic stresses $`\sigma _{ii}=div\sigma `$. Here $`D`$ is the diffusion coefficient and $`\mathrm{\Omega }`$ is the dilatation volume ($`\mathrm{\Omega }>0`$ for interstitials and $`\mathrm{\Omega }<0`$ for vacancies; $`|\mathrm{\Omega }|a^3`$, where $`a`$ is the lattice constant).

In turn, average elastic stresses $`\sigma _{ii}`$ caused by defects satisfy the following equation of the theory of elasticity:

$$\frac{\sigma _{ii}}{𝐫}=3K^{}\mathrm{\Omega }\frac{n}{𝐫},$$
(3)   
where

$$K^{}=K\left(1+\frac{4\mu }{3K}\right)^1,$$
(4)   
$`K`$ and $`\mu `$ are elastic moduli. Substituting Eq.(3) in Eq.(2), we obtain the following expression for $`𝐣`$:

$$𝐣=D\left(1\frac{\mathrm{\Omega }^2K^{}n}{T}\right)\frac{n}{𝐫}.$$
(5)   
It follows from Eq.(5) that the current caused by elastic stresses is directed opposite to the diffusion component of the current. This result has a simple physical meaning. Let us consider, for example, the system of defects with $`\mathrm{\Omega }>0`$. Such defects (e.g., interstitials) stretch the lattice, i.e. increase the average lattice constant, see Eq.(3). But since interstitials themselves tend to move to the stretched regions of the sample, there is an anomalous current along the gradient of defects concentration. This is what we call the deformation interaction of defects.

The steady-state spatially-uniform solutions of Eqs. (1), (5) appear to be unstable at

$$n>n_c=\frac{T}{\mathrm{\Omega }^2K^{}}.$$
(6)   
Hence, in the case that the average concentration of defects exceeds the critical concentration $`n_c`$, the agglomerations of defects (clusters) are formed. For $`T=300`$K and typical values of $`\mathrm{\Omega }510^{23}`$ cm<sup>3</sup> and $`K^{}10^{12}`$ erg/cm<sup>3</sup> we have $`n_c10^{19}`$ cm<sup>-3</sup>, this value being much smaller than the atomic concentration.

We stress that the results presented above are obtained within the framework of a standard defect model in the context of the linear theory of elasticity, without use of any new assumptions and parameters. Note also that the diffusion instability is a common feature of a number of other physical systems, e.g., superconductors and excitonic insulators .

Our results hold true also for defects with $`\mathrm{\Omega }<0`$, e.g., vacancies. Indeed, a vacancy compresses the lattice, i.e. reduce the average lattice constant, while the compressed domain appears to be attractive for other vacancies. If both defects with $`\mathrm{\Omega }>0`$ and $`\mathrm{\Omega }<0`$ (e.g., interstitials and vacancies respectively) are present in the sample, their effective interaction with each other via the lattice deformation is repulsive . Such a repulsion leads to spatial separation of vacancies and interstitials. Interstitials form the interstitial clusters, while vacancies form the vacancion clusters. The interstitial clusters are much less mobile than a single interstitial. As a result, the recombination is suppressed, and the concentration of defects under the influence of an external source can increase monotonously for a very long period of time without saturation . We stress that our model was formulated for the case of a finite defect concentration, and hence does not imply that, e.g., a single vacancy and a single interstitial have a repulsive interaction.

It should be noted that while in isotropic crystals a characteristic range of deformation interaction is very short, about one lattice constant (and hence such an interaction can be viewed as point-like in the continuum limit ), in anisotropic crystals the range of deformation interaction can be much longer. In particular, interlayer interstitials in graphite (both impurities and carbon atoms) are known to deform atomic layers very strongly . As a corollary, the deformation interaction between the defects in graphite should play an important role in temporal evolution of the system of defects under irradiation and influence the spatial arrangement of defects over the sample. The formation of spatially separated interstitial and vacancy clusters and the corresponding suppression of recombination can explain some peculiar features of graphite swelling under irradiation . Note that other mechanisms (e.g., the specific bonding characteristics of the defects) can also contribute to the recombination suppression in graphite.

The effect of interstitial carbon atoms formation in the bulk of graphite has been demonstrated numerically in  making use of empirical interatomic potentials. Below we study the formation of such clusters beneath the surface graphite layer taking into account the strong covalent bonding between carbon atoms constituting the cluster.

3. Computational details

Graphite was modeled by six layers of hexagonal carbon networks, each layer staggered with respect to its neighbors, see Fig.2. The layers were composed of 220 atoms each. The atoms at the periphery of each layer were fixed in order to keep the stability of the crystallite. All other atoms were allowed to move after the interstitials have been inserted into the crystallite. Initially ($`t=0`$) ten interstitial carbon atoms were positioned randomly between one of the peripheral layers (the ”surface layer”) and the penultimate layer.

In the presence of interstitials, the undisturbed graphite layers correspond to a highly nonequilibrium atomic configuration. The potential energy of the system rapidly decreases in time, leading to an increase in the kinetic energy. An additional decrease of the potential energy occurs because of the formation of covalent bonds between the interstitial carbon atoms upon fusion of interstitials into clusters. In order to avoid the overheating of the crystallite, we eliminated the whole kinetic energy of atoms artificially each 70 MD steps (each period of time $`t_o=210^{14}`$ s), i.e. in fact we removed the heat from the crystallite. A rough estimate of the system temperature made on the basis of values of the maximum kinetic energy attained each period $`t_o`$, gives $`T1000`$ K at the very beginning of simulations and $`T10`$ K at the final stage (after the formation of a single interstitial cluster).

We made use of empirical interatomic potentials proposed by Taji et al.  to account for interaction between the carbon atoms located in graphite layers as well as between the interstitial carbon atoms and the atoms of graphite layers. Parameters of those potentials were chosen by fitting the calculated values of interatomic distances and elastic moduli to their experimental values in defect-free graphite. There are two parts in the Taji potential. The first part describes the interaction between the atoms within a graphite layer and takes the non-central force interaction into account in the form of the Keating strain energy potential. This part has an explicit minimum for angles between C-C bonds of 120<sup>o</sup> and hence treats $`sp^2`$ bonding. The second part describes the interaction between graphite layers by the central force potential. This part includes the attractive van der Waals force and the repulsive force that prevents graphite layers from merging together. The second part of the potential was applied by Taji et al. to the interactions not only between the carbon atoms located on the different layers but also between the interstitial carbon atom and other lattice atoms. For an interstitial located between graphite layers, its interaction with graphite layers, as described by the Taji potential, appears to be repulsive, i.e. the interstitial atom does not form covalent bonds with carbon atoms of the nearby graphite layers. Such a choice of the potential agrees with the experimentally observed expansion of graphite due to interstitials that gives evidence for the repulsion between interstitials and graphite layers, i.e. for the absence of covalent bonds between the interstitial atom and the graphite layers. In other words, the interstitial self-energy $`E_{is}`$ is positive. Its value $`E_{is}=`$ 1.8 eV calculated using the Taji potential agrees well with experimental data and theoretical estimates .

Note that the Taji potential was originally proposed  for the case that the interatomic distances and bond angles are not so far from their equilibrium values in graphite. This condition is, however, fulfilled in the simulations of formation of both a single interstitial  and interstitial clusters since the distortion of the surface layer by interstitials is rather smooth, see below.

Of course, the Taji potential is empirical and its use should be also justified by more sophisticated calculations. Since first principles calculations of interstitial characteristics in graphite are not known to us, we have carried out TBMD (tight binding molecular dynamics ) simulations of interstitials between the graphite layers. We have demonstrated that there is indeed a repulsive force between an interstitial and the nearby graphite layers, so that new covalent bonds are not formed and $`sp^3`$ bonding does not appear. We have also shown that there is not only qualitative, but also semi-quantitative agreement between TBMD calculations and simulations based on the Taji potential. In particular, the calculated values of the volume expansion due to an interstitial, the displacements of atoms in the nearest to an interstitial graphite layers, and the value of $`E_{is}`$ are close to those found using the Taji potential.

In order to account for the covalent interaction between the interstitials we employed the TBMD method that had been proven to give a good description of small carbon clusters .

4. Results and discussion

Fig.3 shows the dynamics of 10-interstitial cluster formation in the interlayer region nearest to the graphite surface. Initially ($`t=0`$) the distance between any two interstitials greatly exceeds the characteristic range of covalent interaction ($``$ 2 Å). At the first stage the motion of interstitials is governed entirely by the attractive deformation interaction . At $`t=100`$, a small 2-interstitial cluster is formed. The atoms of the cluster are tightly bound together by covalent bonding, the covalent component of the binding energy being $`E_b=`$ 2.85 eV/atom, in agreement with the experimental value for the dimer C<sub>2</sub> . At $`t=200`$, one more 2-interstitial cluster is formed. In fact, at $`t=200÷500`$ two processes take place, the formation of new diatomic clusters and adsorption of single interstitials by those clusters. The main reason for the motion of single interstitials towards the interstitial clusters is the deformation-induced attraction , as in the case of 2-interstitial clusters formation. At $`t=500`$, there are two 3-interstitial clusters and one 4-interstitial cluster in the interlayer region. Those clusters have the linear chain structure governed by the covalent interaction between the interstitials within the clusters . The covalent component of the binding energy per interstitial is $`E_b=`$ 4.73 eV/atom and $`E_b=`$ 4.75 eV/atom for 3- and 4-interstitial clusters respectively, in accordance with theoretical results of other authors and experimental values (see, e.g., ).

At the next stage of the evolution ($`t=500÷5000`$), three interstitial clusters slowly move towards each other. At $`t=5000`$, two clusters merge into a 6-interstitial cluster. Note that characteristic times of 3- and 4-interstitial clusters formation are an order of magnitude shorter than the time it takes for those clusters to merge together. In other words, the mobility of an interstitial cluster decreases rapidly as the number of interstitials in the cluster increases.

Finally, the 10-interstitial cluster is formed, see Fig.3h. The shape of this cluster is typical for low-dimensional carbon structures that are characterized by the bond angles 180<sup>o</sup> (carbyne) and 120<sup>o</sup> (graphite layers) and by coordination numbers $`Z=2`$ and $`Z=3`$. The covalent component of the binding energy of interstitials within the cluster is $`E_b=`$ 5.50 eV/atom. Note that the most stable geometry of an isolated 10-atom carbon cluster is a ring structure , its binding energy calculated using the same TBMD method equals to $`E_b=`$ 5.87 eV/atom. Nevertheless, we have verified that the shape and the binding energy of the 10-interstitial cluster shown in Fig.3h remained nearly unchanged upon its removal from the crystallite and subsequent relaxation. Hence this cluster, taken alone (outside the graphite) is metastable and corresponds to a local energy minimum.

The deformation of the surface layer by the 10-interstitial cluster is shown in Fig.4. One can see the dome-like feature on the graphite surface. Its height is about 1.5 Å and lateral dimensions are about $`10÷12`$ Å. We suggest that some of the hillocks seen by STM on irradiated graphite surfaces  may be due to relatively large (composed of $``$ 10 carbon atoms) interstitial clusters formed in the interlayer region between the surface and penultimate graphite layers at the sacrifice of deformation interaction of interstitials. The case considered by us (10 interstitials are created simultaneously in the surface area $``$ 10 Å<sup>2</sup>) obviously corresponds to rather high energy of incident ions. Note, however, that interstitials’ dynamics shown in Fig.3 may be viewed as the final stage of cluster formation from well isolated interstitials initially created far apart from each other by low-energy ions and migrated to the region confined by the crystallite due to deformation attraction and/or thermally activated diffusion.

We have made molecular dynamics simulations of interstitial clusters formation for several other initial distributions of interstitials over the interlayer region nearest to the graphite surface. The results are shown in Figs.5-7. They are similar to those presented above. One can see that in either case studied all interstitials, initially located far apart from each other, merge into a single 10-interstitial cluster. The shape of the cluster is not unique, it obviously depends on the initial distribution of interstitials and hence on the specific sequence of small clusters formation and subsequent fusion into a single large cluster. However, the angles between the C-C bonds in those clusters are equal or close to either 180<sup>o</sup> or 120<sup>o</sup>, while the coordination numbers of interstitials within the clusters are $`Z=2`$ or $`Z=3`$ (with obvious exception of boundary interstitials for which $`Z=1`$), just as in the cluster shown in Fig.3h. Thus the shape of the interstitial cluster is strongly influenced by the covalent interactions between the interstitials within the cluster.

The covalent components of the binding energy of interstitials with each other equal to $`E_b=`$ 5.34 eV/atom, 5.32 eV/atom, and 5.68 eV/atom for the clusters shown in Figs. 5b, 6b, and 7b respectively. Those values of $`E_b`$ are lower than the binding energy $`E_b=`$ 5.87 eV/atom of a stable 10-atom ring. Nevertheless we have found that the clusters preserve their shapes upon their removal from the crystallite, the values of $`E_b`$ being changed insignificantly. Thus all final configurations of interstitials are metastable.

Distortions of the surface layer by 10-interstitial clusters are shown in Figs. 5c, 6c, and 7c. It is clearly seen that the shape of the dome-like feature on the graphite surface reflects the particular shape of interstitial cluster. The height of the hillock varies slightly from one interstitial configuration to another, its typical value being 1.5 Å. The lateral dimensions of the hillock are more sensitive to the shape of the cluster, usually being in the range 10 $`÷`$ 15 Å.

Finally, we have also simulated the deformation of the graphite surface by a single carbon interstitial and by interstitial clusters composed of various numbers of interstitial carbon atoms. We have found that the calculated height of the surface hillock increases monotonously with the number of interstitials in the cluster, starting from about 0.8 Å for a single interstitial. The same holds true for characteristic lateral dimensions of the hillock. We believe that such a variation in size of hillocks may be at least a part of the reason for experimentally observed variety of features on irradiated graphite surfaces. For example, Coratger et al.  have found that irradiation of HOPG graphite with 20 keV <sup>12</sup>C<sup>+</sup> ions resulted in protrusions presented an elongated form. These finding can be explained by the formation of relatively small clusters C<sub>N</sub> beneath the graphite surface since for $`N<10`$ a linear chain structure of carbon clusters ia believed either to be the most stable geometry or to lie very close in energy to a ring structure .

5. Summary and conclusions

Making use of molecular dynamics simulations, we have demonstrated that interstitial carbon atoms formed in the interlayer region between the surface and second graphite layers under ion irradiation attract each other. The physical origin for such an attraction is the strong deformation of graphite layers by interstitials that makes it energetically favorable for interstitials to come closer to each other (and thus to minimize the overall buckling of graphite layers) than to stay far apart.

Attractive deformation interaction of interstitials results in formation of interstitial clusters from initially isolated single interstitials. The specific shape of the cluster depends on the initial distribution of interstitials over the interlayer region and is dictated by the strong covalent interaction between the interstitials. Coordination numbers for the majority of interstitials in the cluster are $`Z=2`$ and $`Z=3`$, while all angles between the C-C bonds in the cluster are close to 180<sup>o</sup> or 120<sup>o</sup>, these values being typical for low-dimensional carbon structures. Various arrangements of interstitials in the clusters correspond to different metastable states of the cluster at a fixed number of carbon atoms in it.

Deformation of the surface layer by interstitial clusters results in the dome-like features on the graphite surface. The typical height and lateral dimensions of the hillocks produced by 10-interstitial clusters are 1.5 Å and 10 $`÷`$ 15 Å respectively. Such features have been observed repeatedly by STM on the graphite surfaces irradiated with noble gas ions. The results presented in this paper provide an explanation to experimental observations. In order to make a closer comparison between the theory and experiment, it would be interesting to calculate STM images from the hillocks formed above interstitial clusters of different shapes in graphite.

Acknowledgments

The work was supported by the International Science and Technology Center (Project 467) and by the Russian State Program ”Integration”.

FIGURE CAPTIONS

Fig.1. Distortion of graphite lattice by two interstitials located between the nearest graphite layers at the in-plane distance $`r=`$ 4.9 Å from each other.

Fig.2. Side view (a) and top view (b) of the crystallite composed of 1320 carbon atoms (six layers 220 atoms each). Closed and open circles show the atoms of odd and even (from the top to the bottom) layers respectively.

Fig.3. Dynamics of 10-interstitial cluster formation in the interlayer region nearest to the graphite surface. Top view. Atoms of graphite layers are not shown. Time $`t`$ is measured in units of $`t_o=210^{14}`$ s. $`t=0`$ (a), $`t=100`$ (b), $`t=200`$ (c), $`t=300`$ (d), $`t=500`$ (e), $`t=2000`$ (f), $`t=5000`$ (g), $`t=20000`$ (h).

Fig.4. Deformation of the surface layer by the 10-interstitial cluster.

Fig.5. Initial (a) and final (b) arrangement of 10 carbon interstitials in the interlayer region between the surface and penultimate layers of graphite. Top view. Atoms of graphite layers are not shown. Deformation of the surface layer by the 10-interstitial cluster (c).

Fig.6. Same as in Fig.4 for other initial distribution of 10 interstitials over the interlayer region.

Fig.7. Same as in Fig.4 for other initial distribution of 10 interstitials over the interlayer region.