Atomistic simulations of graphite etching at realistic time scales

Abstract

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 ∼10²⁰ m⁻² s⁻¹. 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.

---

Introduction

Fundamental hydrogen–graphite/graphene interaction has attracted interest in many research fields including astrophysics,¹ 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.¹⁵ A vast amount of knowledge has been gained in fusion energy research, in which dedicated experiments were performed on tokamaks¹⁶ and ion beam setups^17,18 in combination with theoretical studies.^4,19–21 Three fundamental processes were identified: physical sputtering, ion-enhanced chemical erosion, and near-surface sputtering.¹⁵ 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 state^4,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 nanotubes^23–25 and nanocrystalline diamond.²⁶ In graphene/graphite etching simulations, pure graphite or pre-constructed amorphous carbon (a-C:H) samples (to which graphite samples are found to evolve)¹⁵ 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.³⁸ 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,⁴² hydrogen surface diffusion⁴³ or relaxation phenomena,⁴⁴ 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 Carlo^45–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)⁵⁸ 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 ∼10²⁰ m⁻² s⁻¹). 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 package⁶² and modified Colvars module.⁶³ The interatomic potential used in this work is the Reactive Force field (ReaxFF).⁶⁴ In contrast to the widely used 2^nd generation REBO potential,⁶⁵ 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.,⁶⁸ which was previously used to model hydrocarbon desorption and decomposition on Ni nanoparticles as well as CNT growth.⁴⁴ 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 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 (∼10²⁴ m⁻² s⁻¹).^69,70

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.⁷¹

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/(k_bT) and the angle brackets denote the ensemble average, taken as the average over the whole simulation.⁵⁶ 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 (r_i − r_min)/(r_max − r_min). 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 r_max. We used r_max = 2.2 Å and r_min = 1.5 Å for C–C bonds, which is similar to the set-up previously used for hydrocarbon pyrolysis.⁵⁹ 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 r_max 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.

Results

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 sp²/sp³ 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

Fig. 2 Several stages during erosion of the top layer of the graphite substrate: (1) hydrogenation of the top layer, (2) the first CH₂ fragment (red circle), (3) the first CH₃ 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 CH₂ and CH₃ 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 CH₂ and CH₃ 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, CH₂, and CH₃ 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.⁷⁴ This saturation value can partly be explained by the curvature of the graphene layer, which distorts the sp² and sp³ hybridization states due to mechanical stress.⁵ At valley locations this reduces the binding energy and hence incoming H ions are more easily reflected. Once the surface was sufficiently hydrogenated, CH₂ groups appeared which led to an increased H uptake in the first layer. After ∼600 impacts CH₃ 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 sp³ complexes, C–C bond breaking and eventually formation of hydrocarbon complexes at the surface, e.g., CH₃ 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.

Fig. 3 Time evolution of (top) the total hydrogen uptake in the system and the present hydrogen in the CH, CH₂ and CH₃ 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 (C₂H₂ 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.³⁸ Since this repulsion occurs very quickly, the surrounding carbon network has no time to relax.

Fig. 4 Observed erosion mechanisms: (a) hydrogen ion impact induced erosion (the incoming H ion is encircled in window (1), and (b) CH₃ erosion due to a thermal fluctuation.

In the second case, a loosely bound CH₃ 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 sp^x (with free/dangling bonds) which is either caused by the hydrogenation of an sp² C radical or due to abstraction of a bound H atom on an sp³ C atom by an Eley–Rideal type process. At elevated temperature (>400 K), a C atom that contains a hydrocarbon group (e.g., CH₃) and neighbors such an sp^x 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 sp^x C radical (also called β-scission). Alternatively, the hydrocarbon group can directly thermally desorb after the hydrogen irradiation is stopped.²² 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.

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 (C_x>3H_y) 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 10²⁰ m⁻² s⁻¹. 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, CH₂/CH₃ 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 CH₃ 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 CH_x and C₂H_x contributions appear to increase with increasing inter-impact time, while the contribution of large hydrocarbon molecules (C_x>3H_y) 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., C_x=1–3H_y. With regard to the smaller hydrocarbon molecules, the contribution of CH_x is initially lower than that of C₂H_x, 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.

Discussion

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 CH₂/CH₃ group formation is promoted. Additionally, the desorption of small hydrocarbon molecules (C_x=1–3H_y) 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 C₂H_x to CH_x 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 C_x>4H_x 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.³ In this model the total chemical sputtering yield is expressed as:

Y_tot = Y_phys + Y_therm(1 + DY_dam) + Y_surf, (1)

where Y_phys, Y_therm, and Y_surf are the erosion yields due to physical sputtering, chemical erosion and near-surface processes, respectively, and DY_dam is an additional multiplicative term that includes a radiation damage yield Y_dam 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 (T_s) 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 × 10²³ m⁻² s⁻¹ 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.

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 (Y_therm), near-surface erosion (Y_surf), and sp³ concentration (dashed curve, right-axis as indicated by the arrow) are plotted separately. Lastly, the erosion yield and sp³ 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 ∼ 10¹⁹ m⁻² s⁻¹), this regime is already entered for T_s < 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 10²² m⁻² s⁻¹ (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.

Conclusion

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.

Acknowledgements

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.
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.
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.
4. B. V. Mech, A. A. Haasz and J. W. Davis, J. Appl. Phys., 1998, 84, 1655–1669.
5. V. Tozzini and V. Pellegrini, Phys. Chem. Chem. Phys., 2013, 15, 80–89.
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.
7. G. Garberoglio, N. M. Pugno and S. Taioli, J. Phys. Chem. C, 2015, 119, 1980–1987.
8. R. Kumar, V. M. Suresh, T. K. Maji and C. N. R. Rao, Chem. Commun., 2014, 50, 2015–2017.
9. Y. Awano, in Technical Digest – International Electron Devices Meeting, IEDM, 2009.
10. X. Hong, A. Posadas, K. Zou, C. H. Ahn and J. Zhu, Phys. Rev. Lett., 2009, 102, 136808.
11. A. N. Pal and A. Ghosh, Appl. Phys. Lett., 2009, 95, 82105.
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.
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.
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.
15. W. Jacob and J. Roth, Top. Appl. Phys., 2007, 110, 329–400.
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.
17. B. V. Mech, A. A. Haasz and J. W. Davis, J. Nucl. Mater., 1997, 241–243, 1147–1151.
18. B. V. Mech, A. A. Haasz and J. W. Davis, J. Nucl. Mater., 1998, 255, 153–164.
19. J. Küppers, Surf. Sci. Rep., 1995, 22, 249–321.
20. J. Roth and C. Garcia-Rosales, Nucl. Fusion, 1996, 36, 1647–1659.
21. A. Horn, A. Schenk, J. Biener, B. Winter, C. Lutterloh, M. Wittmann and J. Küppers, Chem. Phys. Lett., 1994, 231, 193–198.
22. J. Roth and C. Garcia-Rosales, Nucl. Fusion, 1996, 36, 1647–1659.
23. U. Khalilov, A. Bogaerts, B. Xu, T. Kato, T. Kaneko and E. C. Neyts, Nanoscale, 2017, 9, 1653–1661.
24. U. Khalilov, A. Bogaerts and E. C. Neyts, Carbon, 2017, 118, 452–457.
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.
26. M. Eckert, E. Neyts and A. Bogaerts, Chem. Vap. Deposition, 2008, 14, 213–223.
27. E. Salonen, K. Nordlund, J. Keinonen and C. H. Wu, J. Nucl. Mater., 2001, 290–293, 144–147.
28. E. Despiau-Pujo, A. Davydova, G. Cunge and D. B. Graves, Plasma Chem. Plasma Process., 2016, 36, 213–229.
29. S. J. Stuart, M. Fallet, P. S. Krstic and C. O. Reinhold, J. Phys.: Conf. Ser., 2009, 194, 12059.
30. A. Ito, Y. Wang, S. Irle, K. Morokuma and H. Nakamura, J. Nucl. Mater., 2009, 390–391, 183–187.
31. B. N. Jariwala, C. V. Ciobanu and S. Agarwal, J. Appl. Phys., 2009, 106, 73305.
32. J. Marian, L. A. Zepeda-Ruiz, G. H. Gilmer, E. M. Bringa and T. Rognlien, Phys. Scr., 2006, T124, 65–69.
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.
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.
35. E. D. de Rooij, U. von Toussaint, A. W. Kleyn and W. J. Goedheer, Phys. Chem. Chem. Phys., 2009, 11, 9823–9830.
36. H. Nakamura and A. Ito, Mol. Simul., 2007, 33, 121–126.
37. M. Fallet and S. J. Stuart, Nucl. Instruments Methods Phys. Res. Sect. B Beam Interact. with Mater. Atoms, 2011, 269, 1271–1275.
38. E. Salonen, K. Nordlund, J. Keinonen and C. Wu, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 63, 1–14.
39. E. Despiau-Pujo, A. Davydova, G. Cunge, L. Delfour, L. Magaud and D. B. Graves, J. Appl. Phys., 2013, 113, 114302.
40. P. S. Krstić, S. J. Stuart and C. O. Reinhold, AIP Conf. Proc., 2006, 876, 201–208.
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.
42. Y. Xia, Z. Li and H. J. Kreuzer, Surf. Sci., 2011, 605, L70–L73.
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.
44. E. C. Neyts, A. C. T. Van Duin and A. Bogaerts, J. Am. Chem. Soc., 2011, 133, 17225–17231.
45. E. C. Neyts, B. J. Thijsse, M. J. Mees, K. M. Bal and G. Pourtois, J. Chem. Theory Comput., 2012, 8, 1865–1869.
46. E. C. Neyts and A. Bogaerts, Theor. Chem. Acc., 2013, 132, 1–12.
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.
48. A. F. Voter, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 34, 6819–6829.
49. L. K. Béland, Y. N. Osetsky, R. E. Stoller and H. Xu, Comput. Mater. Sci., 2015, 100, 124–134.
50. E. K. Peter and J.-E. Shea, Phys. Chem. Chem. Phys., 2014, 16, 6430–6440.
51. M. Timonova, J. Groenewegen and B. J. Thijsse, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 144107.
52. H. Xu, Y. N. Osetsky and R. E. Stoller, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 132103.
53. K. M. Bal and E. C. Neyts, J. Chem. Phys., 2014, 141, 204104.
54. M. R. Sørensen and A. F. Voter, J. Chem. Phys., 2000, 112, 9599–9606.
55. D. Perez, B. P. Uberuaga, Y. Shim, J. G. Amar and A. F. Voter, Annu. Rep. Comput. Chem., 2009, 5, 79–98.
56. A. F. Voter, J. Chem. Phys., 1997, 106, 4665–4677.
57. A. F. Voter, Phys. Rev. Lett., 1997, 78, 3908.
58. K. M. Bal and E. C. Neyts, J. Chem. Theory Comput., 2015, 11, 4545–4554.
59. K. M. Bal and E. C. Neyts, Chem. Sci., 2016, 7, 5280–5286.
60. E. C. Neyts and K. M. Bal, Plasma Processes Polym., 2017, 14, 1600158.
61. C. D. Fu, L. F. L. Oliveira and J. Pfaendtner, J. Chem. Theory Comput., 2017, 13, 968–973.
62. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
63. G. Fiorin and J. Hénin, Biophys. J., 2015, 108, 160a.
64. A. C. T. van Duin, S. Dasgupta, F. Lorant and W. A. Goddard, J. Phys. Chem. A, 2001, 105, 9396–9409.
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.
66. A. Ito and H. Nakamura, Commun. Comput. Phys., 2008, 4, 592–610.
67. S. J. Stuart, A. B. Tutein and J. A. Harrison, J. Chem. Phys., 2000, 112, 6472–6486.
68. J. E. Mueller, A. C. T. Van Duin and W. A. Goddard, J. Phys. Chem. C, 2010, 114, 5675–5685.
69. K. Bystrov, L. van der Vegt, G. De Temmerman, C. Arnas and L. Marot, J. Vac. Sci. Technol., A, 2013, 31, 11303.
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.
71. G. J. Martyna, M. L. Klein and M. Tuckerman, J. Chem. Phys., 1992, 97, 2635–2643.
72. S. Casolo, G. F. Tantardini and R. Martinazzo, J. Phys. Chem. A, 2016, 120, 5032–5040.
73. Ž. Šljivančanin, E. Rauls, L. Hornekr, W. Xu, F. Besenbacher and B. Hammer, J. Chem. Phys., 2009, 131, 84706.
74. J. Robertson, Mater. Sci. Eng., R, 2002, 37, 129–281.
75. E. Salonen, K. Nordlund, J. Keinonen and C. H. Wu, Contrib. Plasma Phys., 2002, 42, 458–463.
76. P. Träskelin, K. Nordlund and J. Keinonen, J. Nucl. Mater., 2006, 357, 1–8.

---

Footnote

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

---