Open Access Article
This Open Access Article is licensed under a
Creative Commons Attribution 3.0 Unported Licence

Atomistic simulations of graphite etching at realistic time scales

D. U. B. Aussems *a, K. M. Bal b, T. W. Morgan a, M. C. M. van de Sanden ac and E. C. Neyts b
aDIFFER – Dutch Institute for Fundamental Energy Research, De Zaale 20, 5612 AJ Eindhoven, The Netherlands. E-mail:
bUniversity of Antwerp, Department of Chemistry, PLASMANT Research Group, Universiteitsplein 1, 2610 Antwerp, Belgium
cEindhoven University of Technology, PO Box 513, 5600 MB Eindhoven, The Netherlands

Received 21st June 2017 , Accepted 23rd August 2017

First published on 24th August 2017

Hydrogen–graphite interactions are relevant to a wide variety of applications, ranging from astrophysics to fusion devices and nano-electronics. In order to shed light on these interactions, atomistic simulation using Molecular Dynamics (MD) has been shown to be an invaluable tool. It suffers, however, from severe time-scale limitations. In this work we apply the recently developed Collective Variable-Driven Hyperdynamics (CVHD) method to hydrogen etching of graphite for varying inter-impact times up to a realistic value of 1 ms, which corresponds to a flux of ∼1020 m−2 s−1. The results show that the erosion yield, hydrogen surface coverage and species distribution are significantly affected by the time between impacts. This can be explained by the higher probability of C–C bond breaking due to the prolonged exposure to thermal stress and the subsequent transition from ion- to thermal-induced etching. This latter regime of thermal-induced etching – chemical erosion – is here accessed for the first time using atomistic simulations. In conclusion, this study demonstrates that accounting for long time-scales significantly affects ion bombardment simulations and should not be neglected in a wide range of conditions, in contrast to what is typically assumed.


Fundamental hydrogen–graphite/graphene interaction has attracted interest in many research fields including astrophysics,1 nuclear fusion,2–4 fuel cells,5,6 gas storage,7,8 and nano-electronics.9–14 In particular, research has focused on understanding the release of hydrocarbon molecules from kinetic hydrogen ion bombardment induced chemical reactions, also referred to as chemical sputtering.15 A vast amount of knowledge has been gained in fusion energy research, in which dedicated experiments were performed on tokamaks16 and ion beam setups17,18 in combination with theoretical studies.4,19–21 Three fundamental processes were identified: physical sputtering, ion-enhanced chemical erosion, and near-surface sputtering.15 Using the experimental data and theoretical models as input, the correlation of these processes as functions of quantities such as the ion energy, ion flux, isotope mass, surface temperature and surface state4,16 was described in the semi-empirical Roth–Garcia-Rosales model.21,22

It remains very challenging, however, to confirm the underlying mechanisms of the above mentioned processes on the nano/micro-scale. In this regard, Molecular Dynamics (MD) simulations have been an invaluable tool, not only for the study of graphite/graphene systems as described below, but also for closely related materials such as carbon nanotubes23–25 and nanocrystalline diamond.26 In graphene/graphite etching simulations, pure graphite or pre-constructed amorphous carbon (a-C:H) samples (to which graphite samples are found to evolve)15 are bombarded with energetic hydrogen, thus providing insight into the elementary reactions and emission processes at the atomic level.27–37 For instance, Salonen et al.38 investigated low energy H ion bombardment of a-C:H samples at 300 K, and identified a new sputtering mechanism termed swift chemical sputtering, which could provide a microscopic description of near-surface sputtering. Despiau-Pujo et al.28,39 investigated sputtering under similar conditions, but used multi-layered graphene samples and showed that, due to the lattice structure, surface reactions and erosion become more complex and that before sputtering can occur, initial damage by hydrogenation and vacancy creation has to be induced. MD has also been employed to investigate the mechanisms behind the yield and species composition dependence on quantities such as the ion energy, surface temperature and ion flux.38,40,41

Unfortunately, the ion flux simulated in all of the above mentioned work exceeds the flux range of experiments by at least four orders of magnitude, which directly brings us to the general limitation of MD in terms of accessible time-scales. Typically, the time between two impacts of the etchant species on the surface is in the order of a few ps, after which it is assumed that no further events occur. Practically, in terms of chemical sputtering, this restricts MD simulations to very fast processes only, e.g., physical and near-surface sputtering – both being of the order of 10 fs. Longer time-scale processes (of the order of μs to ms) that become important at elevated surface temperatures and low fluxes, e.g., desorption of weakly bound species,42 hydrogen surface diffusion43 or relaxation phenomena,44 cannot be accessed. It is thus impossible to simulate chemical erosion, and MD studies that aim to find quantitative agreement with experiments under conditions in the range where chemical erosion is dominant have to be considered with caution.

Several methods have been proposed to improve the time-scale reach of atomistic simulations, e.g., (kinetic) Monte Carlo45–53 and accelerated MD methods.54–57 From these methods, hyperdynamics is perhaps the most powerful. In hyperdynamics a bias potential (ΔV) is added to the global potential energy surface (PES) of the system to fill the energy minima, which reduces the waiting time between minima-to-minima transitions. The design of an appropriate bias potential for each system is, however, highly non-trivial. In the recently developed collective variable-driven hyperdynamics implementation (CVHD)58 the bias potential is constructed on the fly in a “self-learning” fashion, allowing the method to be more generically applicable to different systems while requiring little optimization. The flexibility of CVHD is illustrated by the wide range of processes that have already been studied with the method, including surface diffusion, conformational sampling, pyrolysis, combustion, and heterogeneous catalysis, demonstrating its ability to model complex reactions with vastly different reaction rates.58–61

In this work, this method (CVHD) is employed to simulate the erosion of graphite by hydrogen plasma exposure using more realistic inter-impact times (i.e., up to ∼1 ms for an ion flux of ∼1020 m−2 s−1). As a reference, the graphite erosion is first simulated without using a bias potential. Then, the CVHD approach is adopted and the inter-impact time is varied over 9 orders of magnitude. As we will show, this has a significant effect on the ion-induced surface modification and resulting etching rate.

Simulation model

All simulations were performed using the LAMMPS package62 and modified Colvars module.63 The interatomic potential used in this work is the Reactive Force field (ReaxFF).64 In contrast to the widely used 2nd generation REBO potential,65 the ReaxFF potential also includes long-range van der Waals (although this was also added in ref. 28, 66 and 67) and other terms such as Coulomb and torsional interactions. In general, ReaxFF parameterizations are in good agreement with DFT results. In this work we use the parameter set developed by Mueller et al.,68 which was previously used to model hydrocarbon desorption and decomposition on Ni nanoparticles as well as CNT growth.44 This validates its applicability for condensed carbon nanostructures.

The chemical erosion of graphite by the plasma interaction was simulated by impacting the surface with 5 eV H atoms image file: c7sc02763j-t1.tif at random (x, y) positions at normal incidence. The graphite substrate is composed of 4 graphene layers in ABAB stacking, each containing 128 carbon atoms. A sample of a simulation (after 300 H impacts) is depicted in Fig. 1. Periodic boundary conditions were imposed in the x and y directions to mimic a semi-infinite surface. Employing the Nosé–Hoover thermostat, the surface was brought to a temperature of 1000 K, which is the optimal erosion temperature of the maximum ion flux currently achieved in experiments (∼1024 m−2 s−1).69,70

image file: c7sc02763j-f1.tif
Fig. 1 Simulation box of the erosion experiment. Hydrogen is injected at a random (x, y) position with 5 eV at z = 50 Å towards the graphite surface. The lower graphene layer is fixed. Eroded species are removed from the simulation once they enter the “removal zone”. The “freezezone” is introduced to prevent particle loss; in this region particles are not integrated, i.e., the velocity and forces are effectively zero.

After each impact, the motion of each atom was followed for 1 ps in the microcanonical (NVE) ensemble to capture the physics of the hydrogen–surface interaction (e.g., reflection, adsorption or penetration). The timestep throughout the simulation was set to 0.1 fs, which was sufficiently small to conserve the total energy of the system. The integration scheme is the velocity-Verlet algorithm. After the impact of the H ion, the kinetic energy was dissipated into the material (with a typical time constant of 0.1 ps), which caused the average temperature of the substrate to rise. The natural heat conduction out of the cell is mimicked by including an additional canonical ensemble (NVT) phase for a duration of 1 ps. The substrate is cooled to its original temperature by applying a Nosé–Hoover style thermostat with a relaxation constant of 100 fs on all atoms.71

In previous etching simulations that were conducted using MD it would now be assumed that the system remains unperturbed during the residual time before the next impact (∼1 μs to ms) and nothing happens. In this work, however, the CVHD procedure was employed in order to reach the desired longer time-scales. In hyperdynamics, the simulated physical time, also referred to as hypertime, is obtained by multiplying the elapsed MD time by the boost factor 〈eβΔV〉, in which β = 1/(kbT) and the angle brackets denote the ensemble average, taken as the average over the whole simulation.56 The CVHD bias was generated by periodically adding a small repulsive Gaussian potential to the local potential energy landscape at the current state of the system. Thus the gradually increasing CVHD bias is a function of a single collective variable (CV), which includes all degrees of freedom relevant for the process(es) to be observed. A typical example of a collective variable, which was also used in this work, is based on the distortion of bond lengths from their equilibrium values (rirmin)/(rmaxrmin). Care is taken to not add CVHD bias to the transition region, i.e., the region close to the saddle point, as this would corrupt the correct sequence of events in the system, achieved here by choosing an appropriate (i.e., not too large) value of rmax. We used rmax = 2.2 Å and rmin = 1.5 Å for C–C bonds, which is similar to the set-up previously used for hydrocarbon pyrolysis.59 A combination of dynamic and static biasing was used, in which the static base level of the bias potential was set to 0.65 eV, which was found to be below the threshold for an event. On top of this static bias potential, additional Gaussian potentials of width 0.04 and height 0.01 eV were added to the PES every 100 fs, until one of the C–C bonds exceeded rmax for longer than 0.1 ps, after which the bias was removed and a new bias addition was initiated. More details about CVHD and the employed biasing method can be found in ref. 58 and 59. The C–H bond potential should in principle also be biased, but this led to a significant increase in the computation time. Bearing in mind the aim of the current work – to qualitatively show the effect of long time-scales on the physics and chemistry of the simulated system – C–H bond potential biasing will be left for future work.

In contrast to most standard hyperdynamics simulations, no “fixed” time-scale elongation is obtained. That is, the total magnitude of the bias potentials (and hence the overall time-scale that can be reached) is not a fixed quantity, but dependent on the requirements of the system. In addition, each CVHD cycle only lasts as long as needed to reach a certain inter-impact time. Owing to this flexible biasing strategy, arbitrarily long time-scales can be simulated, allowing us to model different ion flux regimes over several orders of magnitude. The typical calculation performance for the simulation with the longest inter-impact interval of 1 ms is ∼1 h wall time (for parallel operation on 4 CPU cores of an Intel® Core i7-2600K, 3.4 GHz, 8 MB cache). From this wall time only up to 66% was spent on CVHD. After each CVHD phase a new impact was initiated (i.e., starting from the NVE phase).

Once erosion was initiated and eroded species entered the removal zone, they were deleted from the simulation after each simulation phase (i.e., NVE, NVT, and CVHD). An additional zone was applied to prevent species from leaving the simulation box; in this zone the atoms were not integrated.


Short time-scale simulation

The graphite erosion was first simulated with an inter-impact time of only 3 ps, from which 1 ps CVHD phase. Fig. 2 shows the time evolution of the top layer of the graphite sample after an increasing number of impacts. From the impacting H ions, a fraction reflected from the surface (mainly from the hollow and bridge sites). Another fraction was chemisorbed on the surface at C–C dimer locations after sp2/sp3 re-hybridization (Fig. 2; 1). This is a similar behavior to that in ref. 28 and is expected because the interatomic potential curves are similarly shaped (for a DFT comparison of the interatomic potential curves near pristine graphene, see Fig. S1). In contrast to ref. 28, however, some H atoms were able to penetrate the top layer and become adsorbed on its back side, e.g., through the bridge site, which may be related to the reduced energy barrier. The chemisorption of the H atoms displays a preference for ortho- and para-pairs – also referred to as clustering, which can be explained by the increased binding energy of these pairs as compared to two separately bound H atoms.39,72,73
image file: c7sc02763j-f2.tif
Fig. 2 Several stages during erosion of the top layer of the graphite substrate: (1) hydrogenation of the top layer, (2) the first CH2 fragment (red circle), (3) the first CH3 fragment (red circle), (4) holes and (5) carbon islands.

The hydrogenation process eventually led to bond breaking due to mechanical stress, either induced by local impacting H ions or by repulsion of a hydrogen pair in the ortho configuration (encircled in Fig. 2; 1). The two unsaturated C atoms that were formed after the C–C bond breaking rapidly saturated their dangling bonds by forming CH2 and CH3 groups, as depicted in Fig. 2; 2–3 (for a DFT comparison of the interatomic potential curves near a defected graphene see Fig. S3). Subsequently, additional H atoms emerged at the back side of the top layer because of surface restructuring. Once the CH2 and CH3 groups were formed, carbon was etched by volatile product formation, as further explained below. With increasing fluence, the C–C rupturing and etching led to holes in the graphene layer (Fig. 2; 4) and eventually to the formation of islands of hydrogenated carbon atoms (Fig. 2; 5), which were etched rapidly. Once holes – so-called etch pits – were formed, the 5 eV H ions could undisturbedly penetrate the active layer and start hydrogenating the second layer. This vertical ‘graphite peeling’ process continued layer by layer, consistent with ref. 28, 30 and 37.

The etching process was investigated more systematically by monitoring the hydrocarbon groups on the surface. In Fig. 3a the hydrogen uptake in the system is shown as a function of the number of impacts as well as the contribution of hydrogen in the CH, CH2, and CH3 groups. It shows that the surface was quickly hydrogenated within 500 impacts and then saturated to a CH/C ratio of ∼40% in the first layer, which is in the range of hard a-C:H films.74 This saturation value can partly be explained by the curvature of the graphene layer, which distorts the sp2 and sp3 hybridization states due to mechanical stress.5 At valley locations this reduces the binding energy and hence incoming H ions are more easily reflected. Once the surface was sufficiently hydrogenated, CH2 groups appeared which led to an increased H uptake in the first layer. After ∼600 impacts CH3 groups were also observed and finally after ∼800 impacts, etching was initiated. This is in line with the chemical erosion model of ref. 19, 21 and 22, in which hydrogenation leads to sp3 complexes, C–C bond breaking and eventually formation of hydrocarbon complexes at the surface, e.g., CH3 radicals. Finally, the loosely bound hydrocarbon groups were etched, as further explained below. Etching continued until all of the carbon and hydrogen atoms in the first layer were released (at ∼2000 impacts). Hereafter, the same cycle was re-initiated.

image file: c7sc02763j-f3.tif
Fig. 3 Time evolution of (top) the total hydrogen uptake in the system and the present hydrogen in the CH, CH2 and CH3 groups, and (bottom) the eroded number of carbon atoms as a function of the number of impacts. The inter-impact time interval is (a) 3 ps and (b) 1 ms.

The etching mechanisms were examined in more detail by observing the exact moment of etch product release. Two possible paths were identified as depicted in Fig. 4. In the first case, an ion impacted very close to a carbon chain (Fig. 4a), eventually causing the release of a hydrocarbon molecule (C2H2 in this case). The release of a weakly bound hydrocarbon molecule due to an ion impact is hereafter referred to as ion-induced erosion. This process points towards swift chemical sputtering as described in ref. 38. In this case, the ions can directly break the covalent C–C bonds of surface hydrocarbon groups bound to the carbon network, because the carbon atoms are forced apart due to the repulsive part of the potential energy function.38 Since this repulsion occurs very quickly, the surrounding carbon network has no time to relax.

image file: c7sc02763j-f4.tif
Fig. 4 Observed erosion mechanisms: (a) hydrogen ion impact induced erosion (the incoming H ion is encircled in window (1), and (b) CH3 erosion due to a thermal fluctuation.

In the second case, a loosely bound CH3 group was desorbed in a thermal fluctuation (Fig. 4b). The release of a weakly bound hydrocarbon molecule due to a thermal fluctuation is hereafter referred to as thermal-induced erosion. The observed mechanism is similar to the ion-enhanced chemical erosion process that was extensively studied in ref. 19, 21 and 22. In these reports, the erosion mechanism is thought to proceed by the following steps. After significant hydrogenation of the surface, a fraction of the carbon atoms will be in an intermediate radical state spx (with free/dangling bonds) which is either caused by the hydrogenation of an sp2 C radical or due to abstraction of a bound H atom on an sp3 C atom by an Eley–Rideal type process. At elevated temperature (>400 K), a C atom that contains a hydrocarbon group (e.g., CH3) and neighbors such an spx C radical can release this hydrocarbon group; the dangling bond of the C atom will form a double bond with the neighboring free bond of the spx C radical (also called β-scission). Alternatively, the hydrocarbon group can directly thermally desorb after the hydrogen irradiation is stopped.22 Remarkably, our simulation shows that this latter process may also occur during hydrogen irradiation, because no hydrogen abstraction was observed before the hydrocarbon complex release.

By determining the exact point at which the hydrocarbon molecule was disconnected from the carbon network, the probability of ion- or thermal-induced release was estimated, as depicted in Fig. 5a. The results show that over 90% of the release is ion-induced.

image file: c7sc02763j-f5.tif
Fig. 5 (a) The probability of the erosion mechanism type (ion- or thermal-induced), for varying inter-impact times: 3 ps, 1 ns, 1 μs and 1 ms. (b) The distribution of the erosion species.

Lastly, the erosion species distribution is depicted in Fig. 5b. Surprisingly, the distribution consists of a significant fraction of large hydrocarbons (Cx>3Hy) compared to the literature.27,30 These large hydrocarbons mainly originate from the carbon islands observed in Fig. 2; 5. Due to the weak van der Waals binding, these groups are only loosely bound to the surface and thus could be desorbed, enhanced by energy transferred from the incoming H ions. These carbon islands predominantly comprise the fraction of hydrocarbons with five or more carbon atoms.

Long time-scale simulation by C–C bond biasing

In the previous section the time between impacts was set to 3 ps. In this section, we apply CVHD to investigate the effect of moving towards more realistic time intervals, e.g., 1 ms, corresponding to an ion flux of 1020 m−2 s−1. In Fig. 3b the time evolution of hydrogen uptake and erosion is plotted for 1 ms inter-impact time intervals. It appears that for all of the time intervals, the erosion process can be described by the same steps as those described in the previous section: hydrogenation, CH2/CH3 group formation, and etching of the top graphene layer. The maximum total H surface coverage – defined as the ratio of the total number of H atoms to the number of C atoms in the top graphene layer – appears to decrease with increasing inter-impact time (from ∼110% to ∼70% H/C) while the relative fraction of CH3 groups increases (from ∼20% to ∼30%). More striking, however, is that the erosion stages appear to occur earlier with increasing inter-impact time. The number of impacts before the H uptake reaches 30% is reduced from ∼370 to ∼150, while erosion starts after ∼80 impacts instead of ∼850 as was the case in the short time-scale simulation. This is also consistent with Fig. 6a, in which the erosion yield (number of eroded carbon atoms per incoming ion) is depicted as a function of the inter-impact time (and ion flux). In the studied range, the erosion yield shows a monotonic increase of a factor of 3.3. For inter-impact times of >10 μs the erosion yield appears to increase more slowly. The erosion yield in the case of a 3 ps inter-impact time (0.07 C/H) is a factor of four higher compared to ref. 28, which is probably related to the lower surface temperature (300 K) in that simulation.

The erosion species distribution for varying inter-impact times is depicted in Fig. 5b. The CHx and C2Hx contributions appear to increase with increasing inter-impact time, while the contribution of large hydrocarbon molecules (Cx>3Hy) decreases. Because the majority of the large molecules are released as carbon islands, this suggests that these clusters disintegrate before they can desorb due to the erosion of small hydrocarbon molecules, i.e., Cx=1–3Hy. With regard to the smaller hydrocarbon molecules, the contribution of CHx is initially lower than that of C2Hx, but it starts to dominate the composition for inter-impact times above 1 ns. Moreover, Fig. 5a shows that between 3 ps and 1 ns a transition occurs from ion- to thermal-induced erosion. In the case of 1 ms inter-impact time, ∼95% of the eroded particles are thermal-induced.


The erosion process is affected in several ways by the long time-scales between impacts. While similar steps appear to be followed as in the short time-scale simulation, these steps are observed to occur earlier in time, i.e., after fewer impacts. This can be understood by considering the biasing acceleration method which has been applied. The biasing effectively brings the C–C bond length close to the value before bond breaking occurs, which simulates the local mechanical stress associated with thermal fluctuations in the graphene surface. With increasing inter-impact time, the effect of these fluctuations will be more pronounced, i.e., for the same number of impacts, C–C bonds are more likely to break (note that this only applies however for the C–C bonds that have an already reduced bond energy, e.g., by the presence of H atoms). The subsequent rise in the number of broken C–C bonds leads to more dangling bonds at the surface, and hence to a higher probability for H ion adsorption. Consequently, saturation is reached earlier while CH2/CH3 group formation is promoted. Additionally, the desorption of small hydrocarbon molecules (Cx=1–3Hy) by thermal fluctuations is enhanced. All these effects lead to a boost in the erosion yield. For inter-impact times exceeding 10 μs, however, the yield appears to increase more slowly. This may be explained by the depletion of the number of potential bond breaking events due to surface relaxation.

Concerning the species distribution, a shift is observed from C2Hx to CHx release with increasing inter-impact time, which may be explained by the transition from ion- to thermal-induced erosion. Moreover, in contrast to other MD work,35,38,40,75 a significant fraction of large hydrocarbon molecules Cx>4Hx is observed in the case of short inter-impact times (3 ps and 1 ns), which we attribute to the release of the carbon islands bound by van der Waals forces. The discrepancy with the literature may be explained by the difference in the sample structure. In the simulations of ref. 35, 38, 40 and 75 an amorphous carbon sample was used instead of graphite, thus carbon islands are not likely to form. Nonetheless, the release of carbon islands may have been boosted due to the shape of the interatomic potential selected in this work. This could result in an overestimation of the erosion yield by less than 60% in the case of a 3 ps inter-impact time Δ.

The trend of a rising erosion yield and the shift towards thermal-induced erosion as a function of the increasing inter-impact time can be explained by the semi-empirical Roth–Garcia-Rosales model.3 In this model the total chemical sputtering yield is expressed as:

Ytot = Yphys + Ytherm(1 + DYdam) + Ysurf,(1)
where Yphys, Ytherm, and Ysurf are the erosion yields due to physical sputtering, chemical erosion and near-surface processes, respectively, and DYdam is an additional multiplicative term that includes a radiation damage yield Ydam and material isotope dependent constant D. Fig. 6a shows the erosion yield calculated by this model for 5 eV ions impacting graphite at a surface temperature (Ts) of 1000 K as a function of the ion flux (φ), hereby adopting the equations and parameters of ref. 3. The results are plotted with and without including the empirically obtained flux effect compensation factor (Cφ−0.54). To gain more insight, the contributions of the thermal and surface erosion (without flux effect) are plotted separately. Note that beyond an ion flux of Γmax = 6 × 1023 m−2 s−1 no experimental data are available and the trend is unknown. The results of the current CVHD study are included along with (for comparison) the data of MD studies in the literature, although different sample structures and temperatures were used. The yields of the literature data were scaled to 5 eV using the ion energy dependence of ref. 70. Three zones can be distinguished, dominated by: (1) desorption, (2) chemical erosion and (3) near-surface-chemical sputtering. Apart from the large deviation between the absolute values of the erosion yields obtained by CVHD and the model, they show qualitatively the same trend in terms of the yield. Additionally, Fig. 5a shows that the transition from ion- to thermal-induced erosion is reproduced. The order of magnitude overestimation of the yield by the CVHD simulation compared to the semi-analytical model can be explained by the difference in rate constants and thus by the precise shape of the interatomic potential, but can also be related to factors which are not considered in the current work, such as the microscopic morphology. Furthermore, the results suggest that the empirically obtained trend of rapidly decreasing erosion yield as a function of the flux (∝φ−0.54) is not observed in our atomistic simulations, at least beyond Γmax, i.e., the data fitting range. In line with ref. 35, this suggests that this effect may not be related to the properties inherent to the material, but is most likely caused by external factors, e.g., redeposition of eroded hydrocarbon molecules or processes occurring in the plasma.

image file: c7sc02763j-f6.tif
Fig. 6 (a) The erosion yield of graphite at 1000 K using our CVHD simulation (●). For comparison, the data from previous MD simulations of H ion bombardment on: a-C:H at 1000 K (▽ (ref. 35)), a-C:H at 300 K (□ (ref. 35), ◊ (ref. 27) and ▷ (ref. 76)) and graphite at 300 K (☆ (ref. 39)). Moreover, the calculated erosion yield using the Roth–Garcia-Rosales model is depicted as a function of the ion flux based on the unmodified model. The contributions of the thermal (Ytherm), near-surface erosion (Ysurf), and sp3 concentration (dashed curve, right-axis as indicated by the arrow) are plotted separately. Lastly, the erosion yield and sp3 concentration are depicted in the case of negligible desorption (dotted curves). (b) The transitional surface temperature from ion- to thermal-induced erosion as a function of the flux (black curve), where the marked area is the regime of dominant thermal-induced erosion. The red curve shows the maximum erosion temperature.

The transition from ion- to thermal-induced sputtering is visualized in Fig. 6b as a function of the surface temperature and ion flux (calculated by the aforementioned model). The marked area shows the regime of dominant thermal-induced erosion, which is inaccessible by common MD simulations. In particular for low flux simulations (Γmax ∼ 1019 m−2 s−1), this regime is already entered for Ts < 670 K. Thus, MD simulation studies on atomic processes in this range have to be considered with caution.

In our simulations C–H bond biasing was not included. Possibly, this led to an underestimation of processes such as thermal hydrogen desorption and hydrogen surface diffusion followed by Langmuir–Hinshelwood recombination. Based on the Roth–Garcia-Rosales model it is expected that neglecting such desorption processes only leads to a saturation of the yield for fluxes below 1022 m−2 s−1 (black dotted line in Fig. 6a), instead of a rapid drop (black dashed line). This trend is consistent with our CVHD results. Nevertheless, unforeseen effects may be significant.


This work presents the effect of long time-scales on graphite erosion by hydrogen ion bombardment using a recently developed hyperdynamics implementation. The results show that while the types of graphite erosion process – hydrogenation, vacancy creation and volatile product formation – do not depend on the inter-impact time, a clear reduction of the required fluence was observed with increasing the inter-impact time. Moreover, the increase in inter-impact time resulted in an increased erosion yield, a reduction of the maximum hydrogen surface coverage and a shift towards smaller hydrocarbon species release. This could be explained by the higher probability for C–C bond breaking due to the prolonged exposure to thermal stress and the associated transition from ion- to thermal-induced etching. In fact, this latter process – chemical erosion – could be accessed for the first time by atomistic simulations due to extended time-scales and is supported by semi-empirical modelling. In conclusion, this study demonstrates that long time-scales can have several important effects on ion bombardment simulations and in contrast to what is typically assumed in the literature these effects may not be neglected, especially for low flux ion bombardment simulations.

Conflicts of interest

There are no conflicts to declare.


DIFFER is part of the Netherlands Organisation for Scientific Research (NWO). K. M. B. is funded as a PhD fellow (aspirant) of the FWO-Flanders (Fund for Scientific Research-Flanders), Grant 11V8915N. The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation – Flanders (FWO) and the Flemish Government – department EWI.

Notes and references

  1. V. Sidis, L. Jeloaica, A. G. Borisov, S. A. Deutscher, F. Combes and G. des Forêts, Cambridge Contemp. Astrophys. Ser., 2000, pp. 89–97 Search PubMed.
  2. G. Federici, C. H. Skinner, J. N. Brooks, J. P. Coad, C. Grisolia, A. A. Haasz, A. Hassanein, V. Philipps, C. S. Pitcher, J. Roth, W. R. Wampler and D. G. Whyte, Nucl. Fusion, 2001, 41, 1967 CrossRef.
  3. J. Roth, A. Kirschner, W. Bohmeyer, S. Brezinsek, A. Cambe, E. Casarotto, R. Doerner, E. Gauthier, G. Federici, S. Higashijima, J. Hogan, A. Kallenbach, H. Kubo, J. M. Layet, T. Nakano, V. Philipps, A. Pospieszczyk, R. Preuss, R. Pugno, R. Ruggiéri, B. Schweer, G. Sergienko and M. Stamp, J. Nucl. Mater., 2005, 337–339, 970–974 CrossRef CAS.
  4. B. V. Mech, A. A. Haasz and J. W. Davis, J. Appl. Phys., 1998, 84, 1655–1669 CrossRef CAS.
  5. V. Tozzini and V. Pellegrini, Phys. Chem. Chem. Phys., 2013, 15, 80–89 RSC.
  6. B. H. Kim, S. J. Hong, S. J. Baek, H. Y. Jeong, N. Park, M. Lee, S. W. Lee, M. Park, S. W. Chu, H. S. Shin, J. Lim, J. C. Lee, Y. Jun and Y. W. Park, Sci. Rep., 2012, 2, 1–6 Search PubMed.
  7. G. Garberoglio, N. M. Pugno and S. Taioli, J. Phys. Chem. C, 2015, 119, 1980–1987 CAS.
  8. R. Kumar, V. M. Suresh, T. K. Maji and C. N. R. Rao, Chem. Commun., 2014, 50, 2015–2017 RSC.
  9. Y. Awano, in Technical Digest – International Electron Devices Meeting, IEDM, 2009 Search PubMed.
  10. X. Hong, A. Posadas, K. Zou, C. H. Ahn and J. Zhu, Phys. Rev. Lett., 2009, 102, 136808 CrossRef CAS PubMed.
  11. A. N. Pal and A. Ghosh, Appl. Phys. Lett., 2009, 95, 82105 CrossRef.
  12. C.-J. Shih, A. Vijayaraghavan, R. Krishnan, R. Sharma, J.-H. Han, M.-H. Ham, Z. Jin, S. Lin, G. L. C. Paulus, N. F. Reuel, Q. H. Wang, D. Blankschtein and M. S. Strano, Nat. Nanotechnol., 2011, 6, 439–445 CrossRef CAS PubMed.
  13. D. Haberer, L. Petaccia, M. Farjam, S. Taioli, S. A. Jafari, A. Nefedov, W. Zhang, L. Calliari, G. Scarduelli, B. Dora, D. V. Vyalikh, T. Pichler, C. Wöll, D. Alfè, S. Simonucci, M. S. Dresselhaus, M. Knupfer, B. Büchner and A. Grüneis, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 165433 CrossRef.
  14. D. Haberer, D. V. Vyalikh, S. Taioli, B. Dora, M. Farjam, J. Fink, D. Marchenko, T. Pichler, K. Ziegler, S. Simonucci, M. S. Dresselhaus, M. Knupfer, B. Büchner and A. Grüneis, Nano Lett., 2010, 10, 3360–3366 CrossRef CAS PubMed.
  15. W. Jacob and J. Roth, Top. Appl. Phys., 2007, 110, 329–400 CrossRef CAS.
  16. J. Roth, R. Preuss, W. Bohmeyer, S. Brezinsek, A. Cambe, E. Casarotto, R. Doerner, E. Gauthier, G. Federici, S. Higashijima, J. Hogan, A. Kallenbach, A. Kirschner, H. Kubo, J. M. Layet, T. Nakano, V. Philipps, A. Pospieszczyk, R. Pugno, R. Ruggieri, B. Schweer, G. Sergienko and M. Stamps, Nucl. Fusion, 2004, 44, L21–L25 CrossRef CAS.
  17. B. V. Mech, A. A. Haasz and J. W. Davis, J. Nucl. Mater., 1997, 241–243, 1147–1151 CrossRef CAS.
  18. B. V. Mech, A. A. Haasz and J. W. Davis, J. Nucl. Mater., 1998, 255, 153–164 CrossRef CAS.
  19. J. Küppers, Surf. Sci. Rep., 1995, 22, 249–321 CrossRef.
  20. J. Roth and C. Garcia-Rosales, Nucl. Fusion, 1996, 36, 1647–1659 CrossRef CAS.
  21. A. Horn, A. Schenk, J. Biener, B. Winter, C. Lutterloh, M. Wittmann and J. Küppers, Chem. Phys. Lett., 1994, 231, 193–198 CrossRef CAS.
  22. J. Roth and C. Garcia-Rosales, Nucl. Fusion, 1996, 36, 1647–1659 CrossRef CAS.
  23. U. Khalilov, A. Bogaerts, B. Xu, T. Kato, T. Kaneko and E. C. Neyts, Nanoscale, 2017, 9, 1653–1661 RSC.
  24. U. Khalilov, A. Bogaerts and E. C. Neyts, Carbon, 2017, 118, 452–457 CrossRef CAS.
  25. U. Khalilov, A. Bogaerts, S. Hussain, E. Kovacevic, P. Brault, C. Boulmer-Leborgne and E. C. Neyts, J. Phys. D: Appl. Phys., 2017, 50, 184001 CrossRef.
  26. M. Eckert, E. Neyts and A. Bogaerts, Chem. Vap. Deposition, 2008, 14, 213–223 CrossRef CAS.
  27. E. Salonen, K. Nordlund, J. Keinonen and C. H. Wu, J. Nucl. Mater., 2001, 290–293, 144–147 CrossRef CAS.
  28. E. Despiau-Pujo, A. Davydova, G. Cunge and D. B. Graves, Plasma Chem. Plasma Process., 2016, 36, 213–229 CrossRef CAS.
  29. S. J. Stuart, M. Fallet, P. S. Krstic and C. O. Reinhold, J. Phys.: Conf. Ser., 2009, 194, 12059 CrossRef.
  30. A. Ito, Y. Wang, S. Irle, K. Morokuma and H. Nakamura, J. Nucl. Mater., 2009, 390–391, 183–187 CrossRef CAS.
  31. B. N. Jariwala, C. V. Ciobanu and S. Agarwal, J. Appl. Phys., 2009, 106, 73305 CrossRef.
  32. J. Marian, L. A. Zepeda-Ruiz, G. H. Gilmer, E. M. Bringa and T. Rognlien, Phys. Scr., 2006, T124, 65–69 CrossRef CAS.
  33. P. S. Krstic, C. O. Reinhold and S. J. Stuart, Nucl. Instruments Methods Phys. Res. Sect. B Beam Interact. with Mater. Atoms, 2009, 267, 704–710 CrossRef CAS.
  34. C. O. Reinhold, P. S. Krstic and S. J. Stuart, Nucl. Instruments Methods Phys. Res. Sect. B Beam Interact. with Mater. Atoms, 2007, 258, 274–277 CrossRef CAS.
  35. E. D. de Rooij, U. von Toussaint, A. W. Kleyn and W. J. Goedheer, Phys. Chem. Chem. Phys., 2009, 11, 9823–9830 RSC.
  36. H. Nakamura and A. Ito, Mol. Simul., 2007, 33, 121–126 CrossRef CAS.
  37. M. Fallet and S. J. Stuart, Nucl. Instruments Methods Phys. Res. Sect. B Beam Interact. with Mater. Atoms, 2011, 269, 1271–1275 CrossRef CAS.
  38. E. Salonen, K. Nordlund, J. Keinonen and C. Wu, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 63, 1–14 CrossRef.
  39. E. Despiau-Pujo, A. Davydova, G. Cunge, L. Delfour, L. Magaud and D. B. Graves, J. Appl. Phys., 2013, 113, 114302 CrossRef.
  40. P. S. Krstić, S. J. Stuart and C. O. Reinhold, AIP Conf. Proc., 2006, 876, 201–208 CrossRef.
  41. L. I. Vergara, F. W. Meyer, H. F. Krause, P. Träskelin, K. Nordlund and E. Salonen, J. Nucl. Mater., 2006, 357, 9–18 CrossRef CAS.
  42. Y. Xia, Z. Li and H. J. Kreuzer, Surf. Sci., 2011, 605, L70–L73 CrossRef CAS.
  43. V. A. Borodin, T. T. Vehviläinen, M. G. Ganchenkova and R. M. Nieminen, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 1–15 CrossRef.
  44. E. C. Neyts, A. C. T. Van Duin and A. Bogaerts, J. Am. Chem. Soc., 2011, 133, 17225–17231 CrossRef CAS PubMed.
  45. E. C. Neyts, B. J. Thijsse, M. J. Mees, K. M. Bal and G. Pourtois, J. Chem. Theory Comput., 2012, 8, 1865–1869 CrossRef CAS PubMed.
  46. E. C. Neyts and A. Bogaerts, Theor. Chem. Acc., 2013, 132, 1–12 CrossRef CAS.
  47. M. J. Mees, G. Pourtois, E. C. Neyts, B. J. Thijsse and A. Stesmans, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 134301 CrossRef.
  48. A. F. Voter, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 34, 6819–6829 CrossRef CAS.
  49. L. K. Béland, Y. N. Osetsky, R. E. Stoller and H. Xu, Comput. Mater. Sci., 2015, 100, 124–134 CrossRef.
  50. E. K. Peter and J.-E. Shea, Phys. Chem. Chem. Phys., 2014, 16, 6430–6440 RSC.
  51. M. Timonova, J. Groenewegen and B. J. Thijsse, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 144107 CrossRef.
  52. H. Xu, Y. N. Osetsky and R. E. Stoller, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 132103 CrossRef.
  53. K. M. Bal and E. C. Neyts, J. Chem. Phys., 2014, 141, 204104 CrossRef PubMed.
  54. M. R. Sørensen and A. F. Voter, J. Chem. Phys., 2000, 112, 9599–9606 CrossRef.
  55. D. Perez, B. P. Uberuaga, Y. Shim, J. G. Amar and A. F. Voter, Annu. Rep. Comput. Chem., 2009, 5, 79–98 CAS.
  56. A. F. Voter, J. Chem. Phys., 1997, 106, 4665–4677 CrossRef CAS.
  57. A. F. Voter, Phys. Rev. Lett., 1997, 78, 3908 CrossRef CAS.
  58. K. M. Bal and E. C. Neyts, J. Chem. Theory Comput., 2015, 11, 4545–4554 CrossRef CAS PubMed.
  59. K. M. Bal and E. C. Neyts, Chem. Sci., 2016, 7, 5280–5286 RSC.
  60. E. C. Neyts and K. M. Bal, Plasma Processes Polym., 2017, 14, 1600158 CrossRef.
  61. C. D. Fu, L. F. L. Oliveira and J. Pfaendtner, J. Chem. Theory Comput., 2017, 13, 968–973 CrossRef CAS PubMed.
  62. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19 CrossRef CAS.
  63. G. Fiorin and J. Hénin, Biophys. J., 2015, 108, 160a CrossRef.
  64. A. C. T. van Duin, S. Dasgupta, F. Lorant and W. A. Goddard, J. Phys. Chem. A, 2001, 105, 9396–9409 CrossRef CAS.
  65. D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni and S. B. Sinnott, J. Phys.: Condens. Matter, 2002, 14, 783–802 CrossRef CAS.
  66. A. Ito and H. Nakamura, Commun. Comput. Phys., 2008, 4, 592–610 Search PubMed.
  67. S. J. Stuart, A. B. Tutein and J. A. Harrison, J. Chem. Phys., 2000, 112, 6472–6486 CrossRef CAS.
  68. J. E. Mueller, A. C. T. Van Duin and W. A. Goddard, J. Phys. Chem. C, 2010, 114, 5675–5685 CAS.
  69. K. Bystrov, L. van der Vegt, G. De Temmerman, C. Arnas and L. Marot, J. Vac. Sci. Technol., A, 2013, 31, 11303 Search PubMed.
  70. J. Roth, R. Preuss, W. Bohmeyer, S. Brezinsek, A. Cambe, E. Casarotto, R. Doerner, E. Gauthier, G. Federici, S. Higashijima, J. Hogan, A. Kallenbach, A. Kirschner, H. Kubo, J. M. Layet, T. Nakano, V. Philipps, A. Pospieszczyk, R. Pugno, R. Ruggieri, B. Schweer, G. Sergienko and M. Stamps, Nucl. Fusion, 2004, 44, L21–L25 CrossRef CAS.
  71. G. J. Martyna, M. L. Klein and M. Tuckerman, J. Chem. Phys., 1992, 97, 2635–2643 CrossRef.
  72. S. Casolo, G. F. Tantardini and R. Martinazzo, J. Phys. Chem. A, 2016, 120, 5032–5040 CrossRef CAS PubMed.
  73. Ž. Šljivančanin, E. Rauls, L. Hornekr, W. Xu, F. Besenbacher and B. Hammer, J. Chem. Phys., 2009, 131, 84706 CrossRef PubMed.
  74. J. Robertson, Mater. Sci. Eng., R, 2002, 37, 129–281 CrossRef.
  75. E. Salonen, K. Nordlund, J. Keinonen and C. H. Wu, Contrib. Plasma Phys., 2002, 42, 458–463 CrossRef CAS.
  76. P. Träskelin, K. Nordlund and J. Keinonen, J. Nucl. Mater., 2006, 357, 1–8 CrossRef.


Electronic supplementary information (ESI) available. See DOI: 10.1039/c7sc02763j

This journal is © The Royal Society of Chemistry 2017