Simulating Radiation Damage Cascades in Graphite

A B S T R A C T Molecular dynamics simulation is used to study radiation damage cascades in graphite. High statistical precision is obtained by sampling a wide energy range (100–2500 eV) and a large number of initial directions of the primary knock-on atom. Chemical bonding is described using the Environment Dependent Interaction Potential for carbon. Graphite is found to exhibit a radiation response distinct from metals and oxides primarily due to the absence of a thermal spike which results in point defects and disconnected regions of damage. Other unique attributes include exceedingly short cascade lifetimes and frac-tal-like atomic trajectories. Unusually for a solid, the binary collision approximation is useful across a wide energy range, and as a consequence residual damage is consistent with the Kinchin–Pease model. The simulations are in agreement with known experimental data and help to clarify substantial uncertainty in the literature regarding the extent of the cascade and the associated damage. Despite being one of the original nuclear materials, remarkably few molecular dynamics simulations have been performed to understand radiation response in graphite. Whereas a vast computational literature exists for radiation processes in metals and oxides (see Refs. [1–4] for reviews) only a handful of simulations exist for graphite due to historical difficulties associated with describing bonding in carbon [5]. Aside from point defect energetics [6,7], and threshold displacement energies [8–10], little is known atomistically about cascade behavior, recovery following ballistic displacement or temperature-driven dynamical effects. In a modern context, understanding of radiation processes in graphite is motivated by lifetime extensions of existing graphite-moderated reactors [11,12], and future Generation-IV technologies such as the high-temperature graphite-moderated design [13]. The simulation of radiation damage using molecular dynamics (MD) has a long history, extending back to the first ever MD publication in 1960 on focussed collision sequences in copper [14]. A great number of radiation cascade simulations were performed over the following decades, facilitated by the development of the embedded atom method [15,16] for metals and Buckingham-type potentials [17,18] for ionic solids and oxides. At first glance many cascades look the same, commencing with a ballistic phase in which the kinetic energy of the primary knock-on atom (PKA) is rapidly deposited into a small region, followed by transient localised


Introduction
Despite being one of the original nuclear materials, remarkably few molecular dynamics simulations have been performed to understand radiation response in graphite. Whereas a vast computational literature exists for radiation processes in metals and oxides (see Refs. [1][2][3][4] for reviews) only a handful of simulations exist for graphite due to historical difficulties associated with describing bonding in carbon [5]. Aside from point defect energetics [6,7], and threshold displacement energies [8][9][10], little is known atomistically about cascade behavior, recovery following ballistic displacement or temperature-driven dynamical effects. In a modern context, understanding of radiation processes in graphite is motivated by lifetime extensions of existing graphite-moderated reactors [11,12], and future Generation-IV technologies such as the high-temperature graphite-moderated design [13].
The simulation of radiation damage using molecular dynamics (MD) has a long history, extending back to the first ever MD publication in 1960 on focussed collision sequences in copper [14]. A great number of radiation cascade simulations were performed over the following decades, facilitated by the development of the embedded atom method [15,16] for metals and Buckingham-type potentials [17,18] for ionic solids and oxides. At first glance many cascades look the same, commencing with a ballistic phase in which the kinetic energy of the primary knock-on atom (PKA) is rapidly deposited into a small region, followed by transient localised melting in the form of a thermal spike, and concluding with rapid cooling and partial recrystallization. Depending on the material type and chemistry, the entire process results in disordered structures such as stacking fault tetrahedra, isolated point defects, amorphous pockets or extended defect clusters. For highly radiation tolerant materials (e.g. rutile TiO 2 [19]) it is even possible to have no residual defects at all, such is the extent of the driving force towards the crystalline state. Regardless of the outcome, the behavior is typically dictated by an interplay between the native crystal structure and dynamic recombination/annealing of defects.
To the best of our knowledge, the first reported simulation of radiation cascade effects in graphite was performed in 1990 by Smith [20], who used the Tersoff potential [21] to study self-sputtering and related phenomena. Further work by Smith and Beardmore [8] expanded the computational techniques to include potentials proposed by Brenner [22] and Heggie [23] and examined bombardment with Ar and C 60 . Key results included quantification of ion-surface interactions and an estimate of the threshold displacement energy E d of 34.5 eV. Similar studies were performed by Nordlund et al. [24] using a long-range extension of the Tersoff potential to quantify defect creation responsible for hillocks observed at the surface of graphite. Their potential has been used to study ion irradiation in a variety of sp 2 carbon systems (for a comprehensive review see [25]) but it has not been applied to damage cascades in graphite. In work motivated by next-generation reactor design, Hehr et al. [10] modified the Brenner potential to study temperaturedependance of E d , finding values of 44.5 eV at 300 K and 42.0 eV at 1800 K. They did not, however, report any calculations of radiation damage cascades.
Here we report graphite cascade simulations using the Environment Dependent Interaction Potential (EDIP) for carbon [26,27] coupled with the standard Ziegler-Biersack-Littmack (ZBL) potential [28] to describe close-range pair interactions. Originally developed to study thin film deposition of amorphous carbon, EDIP has since been applied to study numerous other carbon forms including carbon onions [29,30], glassy carbon [30,31], peapods [32], nanotubes [30,33] and nanodiamond [34]. Our article is structured as follows: in Section 2 we detail our methodology for linking the EDIP and ZBL potentials and outline our procedure for performing simulations and defect analysis. In Section 3 we consider first the qualitative behavior of cascades in graphite, considering specific examples which illustrate binary-collision-type behavior and channeling. This is followed by quantitative analysis averaged over a large number of PKA directions, examining cascade properties such as timescale, length scale and defect production. We conclude in Section 4 with a discussion of how radiation damage in graphite is fundamentally different to metals and oxides and link our results with historical models in the literature.

Methodology
The radiation damage cascades are simulated using MD with equilibrium interactions governed by the EDIP methodology for carbon [26,27]. The potential is based on an earlier EDIP potential for silicon [35], where the key elements are two-body and three-body interactions modulated by an atomic bond-order term derived from the coordination. The carbon variant of EDIP includes a more sophisticated aspherical coordination counting term which provides an excellent description of bond-making and breaking, in particular the energy barrier for conversion between sp 2 and sp 3 hybridization such as occurs in the graphite/diamond transformation. Its ability to also accurately describe disordered states makes EDIP ideal for studying radiation damage.
To accurately model close-ranged interactions between atoms, the pair potential within the EDIP formalism is smoothly switched to the ZBL pair potential [28] at small separations. Due to the environmental dependence, this transition is less straightforward than with pair potentials where interpolation functions such as cubic splines can be used to connect the potentials. Our approach uses two Fermi-type scaling functions (SF EDIP , SF ZBL ) which are defined using three parameters; the positions of the midpoint of each function (r EDIP and r ZBL ) and the width w of the switching region. After a trial-and-error process, values of r EDIP ¼ 1:05 Å , r ZBL ¼ 0:45 Å , and w = 0.07 Å were selected to ensure a smooth transition between the two regimes and to avoid inflexion points associated with changes in curvature. A plot of the switching functions along with the effective pair potential for a coordination number of three are shown in Fig. 1. Note that this is the same approach as employed in previous EDIP simulations of ion impact [36,37] where the scaling function approach was first introduced and briefly defined.
Simulations were carried out in lattices equilibrated at 300 K and after the PKA was initiated the motion was followed for 5 ps. At the conclusion of the simulation steepest descent minimization was performed prior to defect analysis. Periodic boundary conditions were employed in each of the cartesian directions along with a 3.5Å thermal layer to remove excess kinetic energy during the cascade simulations. In addition, a layer of atoms perpendicular to the basal plane were held fixed to prevent any net transverse motion of the planes. To produce a representative set of collision cascades, a range of PKA energies and directions of initial velocity were sampled. Cascades were initiated with PKA energies of 100, 250, 500, 750, 1000, 1500, 2000 and 2500 eV within orthorhombic supercells as large as 157.7 · 152.6 · 153.6 Å 3 and containing up to 440,448 atoms. As in the work of Robinson et al. [38], the initial directions of the PKA were determined by uniformly distributing points on a unit sphere (the so-called Thomson problem [39]) to select the direction of the initial velocity. This methodology allows for an arbitrary number of directions to be chosen. The data presented here is a composite of two uncorrelated data sets, one with 10 uniformly distributed data points and a second containing 20 points. To assist in retaining the cascade within the simulation cell, the location of the PKA atom was chosen such that the direction of initial motion was towards the center of the supercell. Cascades were analysed using various quantitative measures. The number of defects were computed using a vacancy radius v r of 0.9 Å [40]. This value was also used to determine atomic displacements. Coordination numbers are determined using a nearest-neighbor cutoff of 1.85 Å . At higher PKA energies the graphite planes are prone to a degree of buckling due to the weak interlayer interactions. In a number of simulations this proved problematic, with significant numbers of atoms being identified as defects by the algorithm, even though their local environment was still purely graphitic. To circumvent this problem we added an additional criteria in which we computed the displacement perpendicular to the basal plane for all 3-fold coordinated atoms. If the magnitude of this displacement was less than 1.8 Å these atoms were identified as crystalline and removed from any subsequent defect analysis. Visual inspection of the problematic data sets showed that this approach worked extremely well, and meant the automatic defect counting algorithm could be used with confidence.

3.1.
Individual cascades Fig. 2 shows a representative example of a 1 keV cascade in graphite. Blue squares denote vacancies while red circles denote interstitials. In the first frame (t = 0.002 ps) the PKA has moved on the order of the vacancy radius and has left behind a vacancy; the initial direction of the PKA is indicated by the arrow. After 0.014 ps the PKA has already experienced a close-approach collision with another carbon atom, effectively splitting into two sub-cascades. In the third and fourth frame (t = 0.03 and 0.04 ps) it is clear that the upper sub-cascade is the more energetic of the two and continues to exhibit a branching structure with each successive collision. In contrast, the less energetic sub-cascade is already close to its end-of-range point. By 0.06 ps the entire cascade is close to maximum extent, with most of the remaining excess kinetic energy concentrated in a handful of atoms. At 0.13 ps the atoms no longer have sufficient energy to create new defects and self-induced annealing becomes the dominant process. The relaxation process is extremely rapid, and the structure at 1.40 ps is essentially identical to that when the simulation concludes. By this stage the number of defects has decreased considerably, from a peak of 37 interstitials at 0.11 ps, to 10 at 1.40 ps.
The final frame at the bottom-right of Fig. 2 shows an alternative view of the final structure, involving a rotation to highlight the cascade relative to the graphitic planes and a color-coding to indicate variations from the standard coordination number of three. At the bottom left a single undercoordinated atom results from the vacancy created by the PKA, while the green overcoordinated atoms mostly arise in the graphitic planes adjacent to an interstitial atom position between two-layers, thereby increasing the coordination number of the in-layer atoms to four. This is particularly apparent for the defect complex at the top of the panel where the yellow trajectory trace follows the path of the mobile atom which has moved between adjacent layers for a short distance, eventually creating an interlayer defect. Density functional theory calculations [41] have identified a closely related configuration, known as the spiro interstitial due to its resemblance to the spiro-pentane molecule. Some of the structural details of the spiro interstitial (specifically, the two sets of triangular C-C bonds) differ to the MD simulation, but the origin of the discrepancy is well-known and arises from neglect of three-center terms, a common approximation in pair potentials and tight-binding methodologies [42]. The key observation is that the interlayer defect is indeed prevalent in the simulations and is strongly correlated with overcoordinated atoms created by interstitials.
One of the striking features in Fig. 2 is the fractal-like branching structure of the defect trajectories. This arises from the highly energetic collisions involving the PKA and subsequently displaced lattice atoms. To convey a sense of the strength of these interactions we computed the maximum kinetic energy of any atom, KE max , and plotted this quantity as a function of time (see Fig. 3). The arrows in the figure highlight the moments of close approach in which the PKA (or another more energetic atom) interacts strongly with an atom on a lattice site. At the instant of closest approach the velocities are transiently very small and the forces are enormous; the repulsion between the two atoms then splits the cascade and converts potential energy into kinetic energy spread between the two atoms. Correspondingly, each branching point seen in the time-series snapshots of Fig. 2 can be correlated to one of the close-approaches denoted by the arrows in Fig. 3.
Consideration of the details in Fig. 3 reveals an aspect which is initially surprising, namely, an extremely short time period when atoms have substantial kinetic energy. Cascade simulations in metals and oxides with PKA energies in the keV range typically evolve on the timescale of a picosecond or so, and yet here the maximum kinetic energy has fallen below 10 eV after just 0.090 ps. To quantify this behavior we counted the number of atoms with a kinetic energy exceeding thresholds of 1 and 10 eV. The latter is a measure of the number of ''fast atoms'' which might be reasonably expected to be associated with motion of atoms to a different lattice site or defect configuration, while the lower threshold conveys a sense of the rate at which the high thermal conductivity of the surrounding matrix removes heat from the cascade. Fig. 4 plots the time-dependence of both quantities and confirms the earlier impression from Fig. 3. The maximum number of fast atoms (red trace) is achieved after just 0.032 ps, while the number of ''warm atoms'' (blue trace) first reduces to zero at 0.235 ps. These two quantities are indicated by arrows in Fig. 4, and in preparation for later use we denote them t max and t end , respectively. This exceptionally rapid time evolution is common to all cascades considered in this study and is quantified statistically in the following section. Given the striking difference compared to other materials, we calculated the speed of sound in a diamond rod, [43] and reproduced the known experimental value of around 12 km/s. This confirms that the rapid dynamics of the cascade is a bonafide effect and highlights the value of studying graphite as a contrasting materials system as compared to the well-known radiation effects in metals and oxides.
On a technical level, the short-lived but highly energetic events in Figs. 3 and 4 highlight the importance of the variable timestep algorithm. When the PKA is initiated the equations of motion are well-integrated using a timestep of 0.018 fs, but during the first collision this falls to just 0.0022 fs for a brief period, before increasing to a maximum of 0.025 fs midway between the first and second collisions. With each collision the timestep is temporarily reduced and the cycle repeats, and the timestep gradually trends upwards towards a constant value of 0.23 fs by the end of the simulation. Across each collision the conservation of energy is excellent, leading to shifts no greater than 0.1 eV and typically far less. Due to the efficiency of the variable timestep algorithm, the total simulation of 5 ps was completed in fewer than 30,000 timesteps, with around 10,000 steps required to cover the first picosecond.
A second example of a 1 keV cascade in graphite is shown in Fig. 5(a), this time employing a different initial PKA direction. In this instance the cascade proceeds in a manner very different to that in Fig. 2. A branching structure is not observed, and instead the PKA is deflected into a channeling direction with a 10 12 orientation. Once in the channel the PKA travels a substantial distance without creating permanent defects, losing energy at a constant rate of 18 eV per layer traversed. The kinetics of this process are summarized by Fig. 5(b) which shows the time evolution of the most energetic atom in the cascade. Up until an end-of-range collision at 0.083 ps, the atom in question is always the PKA. Whilst in the channel the PKA undergoes a series of collisions in which it passes through the middle of a hexagonal ring of atoms, losing kinetic energy on entry and regaining a portion of this upon exit. Since the energy loss of 18 eV/layer is below the threshold displacement energy (which we show below to be around 25 eV), the disruption is transient and the lattice quickly recovers. The process is repeated nine times, during which the PKA travels more than 30 Å and loses around 160 eV without creating a single permanent defect.
The channeling direction along which the PKA travels is shown in ball-and-stick form in the inset within Fig. 5(b). It is immediately apparent that this is not a prototypical channel such as found in a lattice with cubic symmetry, but instead a pseudo-channel in which the PKA must follow an undulatory path as it progresses through one layer to the next. We are not aware of this channeling direction having previously been identified for graphite, a situation that perhaps reflects the difficulty of anticipating the pseudo-channel in the first place. Even for conceptually straightforward channels such as those parallel or perpendicular to the c-axis, there are significant experimental difficulties associated with preparation of the sample and its alignment relative to the incident beam [44]. We return to this question of channeling in Section 4 where we discuss the simulations in the context of the historical literature.

Statistical analysis
Having outlined some of the qualitative features of radiation cascades in graphite, we now proceed to a quantitative treatment of various key properties. Statistical sampling is a crucial element of any cascade simulation analysis, and particularly so for graphite where the low packing fraction and anisotropic crystalline structure facilitates strong variations in radiation response as a function and direction. Up to thirty directions are sampled for each PKA energy, providing a high degree of precision which enables extraction of clear trends in the data. Fig. 6 shows the energy dependence of the quantities t max and t end defined in Fig. 4. The solid points show the average value while the error bars denote one standard deviation. The magnitude of the latter indicates the high degree of variability between individual cascade events and highlights the importance of sampling many uncorrelated directions. Due to the large number of simulations, the standard error in the mean is around a factor of 5 smaller than the ranges shown in the figure, and hence the solid points provide a good estimate of the true mean of both quantities. One of the surprising facts to emerge from Fig. 6 is the weak energy dependence  of the time of peak defects (upper panel), and spike cooling time (lower panel). For example, increasing the PKA energy by a factor of four from 500 eV to 2 keV changes t max by only 70% and t end by 55%. This behavior reflects the fractal-like branching structure seen earlier in Fig. 2 where the cascade behavior is almost entirely ballistic, involving repeated splitting into a series of sub-cascades. Accordingly, increasing the PKA by a large factor makes only a relatively small difference to the cascade lifetime.
To quantify the difference between graphite cascades and those in most other materials we show in the lower panel of Fig. 6 the energy dependence for a thermal spike produced when the energy of the PKA is delivered into a relatively compact region and induces local melting. Presuming a spherical spike for simplicity, an analytic solution [45] of the heat diffusion equation shows that the cooling time varies as E 2=3 , where E is the PKA energy. This energy variation is shown as a solid black line in Fig. 6, using the data point at 500 eV as an arbitrary anchor point. Clearly the graphite cascades have an energy dependence entirely unlike the thermal spike model, confirming the previous visual impressions seen in Figs. 2 and 5 that localised melting and extended disordering do not occur in graphite.
Having obtained these two numerical data sets, we subsequently explored various mathematical functions, and found that both data sets are closely described by a power-law expression of the form aE x , were E is the energy of the PKA. The fitted functions, shown as a blue-dotted line in each panel, reproduce the numerical data extremely well over a wide energy range, extending even to relatively low energies approaching the threshold displacement energy. For t max , which is the noisier of the two data sets due to smaller numerical values, the exponent x is 0.37, while for the cascade lifetime t end the value x is 0.28 and the fitting quality is excellent. Since the exponent is so small, increasing the PKA energy to much larger values, for example, to 20 keV, only increases the spike lifetime by 80% relative to a 2.5 keV cascade. While we do not presently have a physical interpretation for this power-law behavior, the quality of the fit is striking and it would be instructive to test the trend for much higher energies and to examine whether it can be exploited in simplified, non-MD, models of cascade evolution. Fig. 7 shows the energy dependence of the size of the cascade and the range of the PKA. The latter is determined by simply taking the difference between the initial and final positions of the PKA atom, while the cascade length is defined as the largest distance between any two defects in the cascade. For the purposes of this analysis, a defect was a defined as any atom with a local potential energy greater than À7 eV/ atom. Since the cohesive energy of graphite in the EDIP formalism is À7.361 eV/atom, this corresponds to a strain energy of around 0.35 eV/atom. Visual inspection of the atoms identified in this way showed this measure provided a simple and accurate measure of the region affected by the cascade. As in Fig. 6, each data point is an average across as many as 30 different directions and the error bars denote one standard deviation. Also shown in Fig. 7 are straight line fits to the two data sets. To a good approximation both the PKA range and cascade extent scale linearly with energy. The gradient of both quantities are quite similar, 27 Å /keV for the PKA range and 31 Å /keV for the cascade extent. For the highest energy cascade of 2.5 keV the PKA range is around 85% of the cascade extent, and slightly smaller at lower energies. An immediate consequence of this linear variation is that increasing the PKA energy comes at high computational cost; numerical specifics and possible alternative strategies are outlined in the Discussion.
To quantify the process by which the PKA energy is converted into displacements we tracked the maximum kinetic energy KE max in every simulation and analyzed its time dependence. Motivated by data such as that shown in Fig. 3, a collision was recorded if KE max increased by more than 5 eV following a minimal value. This metric robustly identifies all of the heavy collisions indicated by arrows in Fig. 3, as well as the much larger number of collisions for the channeling process in Fig. 5. Once KE max becomes sufficiently small, circa 50-100 eV, collisions are no longer identified since a subsequent rise in kinetic energy does not occur following a close approach.
The results of the analysis are summarized in Fig. 8 error bars indicating one standard deviation. The straight line is a linear fit to the data and intersects the average values with high accuracy; a companion plot employing the standard error of the mean is not shown as the resultant uncertainties are comparable in size to the circles used in the figure. Considering a PKA energy of 1 keV as an example, we see that on average less than 10 collisions are required to thermalize the PKA, and that accordingly the branching cascade seen in Fig. 2, with 9 collisions, is more typical than the channeling cascade (Fig. 5) which thermalizes after 19 collisions. From the inverse slope of the linear fit to the data we see that a typical collision loses 112 eV, a result which is consistent with visual inspection of the branching cascade. The linear behavior holds over the entire energy range considered, and makes a useful starting point for estimating the behavior of more energetic cascades.
Two of the most crucial quantities to quantify in a radiation cascade simulation are the number of displaced atoms and the number of defects created. As noted in Section 2, we define defects and displacements using a vacancy radius of 0.9 Å . The energy dependence of both data sets is shown in Fig. 9, along with the theoretical behavior predicted by the Kinchin-Pease (KP) [46] and Norgett-Robinson-Torrens (NRT) [47] models. Both models employ the threshold displacement energy, E d , as a single parameter; in the KP model the number of displacements is computed as E=ð2E d Þ, where E is the PKA energy, while in the modification of NRT an additional multiplicative factor of 0.8 is included to empirically describe recombination. Consistent with the experimental and computational literature, a value of E d ¼ 25 eV was assumed for the two models. The simulations show excellent correlation with the two theories, a result which is somewhat remarkable in itself given that the theoretical underpinnings of the KP and NRT treatments are not especially strong, especially regarding the NRT factor of 0.8 which dates back to a computer simulations performed in the mid 1970's. In many metal and oxide systems defect recombination in the post-ballistic phase can lead to vastly fewer defects than displacements, and hence the KP/NRT methods don't always provide predictive power. For graphite, however, these simple models work extremely well, even down to the empirical factor of 0.8 which relates displacements to defects.

Discussion and conclusion
One of the main insights to emerge from this study is the profound difference between cascades in graphite as compared to other widely-studied solids. The branching structure, absence of localized melting, and KP-type defect generation are all examples of behavior which place graphite in a special category. Equally interesting is that many of these insights were correctly qualitatively understood long ago, as demonstrated by Fig Both schematics date from the 1960's, with panel (a) from Nightingale [48] and panel (b) from Simmons [49]. Although one could argue the minor details of the schematics, the broad picture is clearly correct, particularly regarding the branching structure of the trajectories. As for the origin of this behavior, a definitive answer cannot yet be given, but possible reasons include the low mass of carbon relative to many solids, the high thermal conductivity of graphite, and the low-packing fraction. All of these ideas are amenable to computer simulation through controlled comparison studies, and represent a promising direction for future work.
Where the simulations offer a clear advance over previous empirical understanding is in quantitation. Key quantities such as the threshold displacement energy, PKA range, mean free path, etc, have been subject to considerable uncertainty. In the case of the range of the PKA, Simmons [49] quotes a value of 67 Å for a 1 keV PKA, as compared to the 30 Å computed here. The same source similarly overestimates energy loss per collision, listing a value of 196 eV per collision as compared to our value of 112 eV per collision. Earlier hardsphere results reported in Nightingale [48] disagree to an even larger extent, predicting a mean-free-path between collisions of 84 Å for a 1 keV PKA. Establishing E d for graphite has also been highly problematic, with literature estimates covering a wide range, spanning 10 to 60 eV [50,51]. In contrast, our value of 25 eV determined by the KP and NRT relationships is statistically sound. Burchell has previously noted [52] that a value of 60 eV has gained wide acceptance (specifically in the nuclear industry) but argues that a much smaller value of 30 eV is appropriate; the estimate here of 25 eV is broadly consistent with this assessment. We note that Yazyev et al. [9] also estimated an E d value of 25 eV from their density functional theory calculations, while Smith and Beardmore [8] reported 34.5 eV, and Hehr et al. [10] reported around 45 eV. As a caveat, we note that direct analysis of MD trajectories (such as in [38] can provide an even more refined estimate of E d , and so the present number should not be considered the final value. Such an approach is the most unambiguous route to determining E d , since the analysis is based entirely on counting statistics over a large number of trajectories, and no inputs or prior functional form need to be assumed, save the interatomic potential itself. Regarding channeling, there have been a few experimental studies [44,53], but the measurements are difficult and sample preparation and alignment are paramount. Channeling down a h0 0 0 1i channel (i.e. parallel to the c-axis) has been observed (see Fig. 1 in [53] for a schematic), and long ago it was proposed [54] that channeling might occur in the h1 1 2 0i direction. The simulations show that while channeling is not a common occurence, it is certainly possible under certain conditions. To the best of our knowledge the channeling process shown in Fig. 5 has not been previously described or envisaged and is an excellent example how computer simulation can provide new insights into radiation phenomena. The closest connection with previous work is the calculation of Kaxiras and Pandey [55] which found that a hexagon center defect in graphite has an energy of 19.5 eV, close to the energy loss rate observed in the simulations, suggesting that the PKA has to climb up the potential hill of $18 eV before being squeezed out again. Finally, we note that one of the factors which limits channeling is that the atoms are invariably initially located on a lattice site and hence displacement tends to cause an immediate collision with a nearby atom. A different situation exists for interstitial carbons, but such atoms are present in low numbers and cannot be considered typical.
Looking to future simulations of cascades in graphite, several conclusions immediately present themselves. Firstly, there are promising prospects to study how temperature affects the evolution of the cascade and the associated relaxation dynamics. As discussed by Kelly [56], a large experimental literature exists for graphite on dimensional change, mechanical behavior and thermal properties. However, simulations have not been applied to the atomistic perspective. It would be particularly fruitful to combine defect analysis as carried out in this work with density functional theory calculations of defect energetics and activation barriers for migration and recombination. This knowledge would be particularly beneficial for understanding graphite structural evolution under irradiation such as the recently proposed ruck-and-tuck model [57]. On a technical level, one unexpected detail is that wall thermostats are not essential as the simulation cells are sufficiently large that the excess kinetic energy of the PKA is easily accommodated within the cell. For the simulations performed here, the excess energy was typically 0.005 eV/atom and hence the temperature rise is of the order tens of Kelvin. More significant is the very large number of atoms required to contain the cascade as the PKA energy rises. Extrapolating the data in Fig. 7 to higher energies shows that very large simulation cells are required to reliably avoid the cascade interacting with the boundaries. To contain a 5 keV cascade at a 2r confidence level (97.5%) requires a cell of 2 million atoms; the required supercell side-length is 263 Å , comprising 253 Å from the extrapolation of the mean cascade length plus two standard deviations, and a boundary layer of 10 Å . Higher energies make the numerics even more extreme, as the number of atoms required scales as E 3 due to the linear variation shown in Fig. 7. For example, 40 keV cascades have been simulated in a variety of oxides [58], but the same calculations in graphite would require 1 billion atoms, far beyond what is presently practicable for carbon. For large systems it may be preferable instead to develop a stochastic approach based upon the atomistically-derived information extracted from molecular dynamics. Such a scheme should in-principle be possible given the high statistical reproducibility evident in the computational data presented here.
In summary, we have performed molecular dynamics simulations of radiation damage cascades in graphite. To the best of our knowledge, this is the first comprehensive study of its kind. We find strikingly different behavior to metals and oxides, with the graphite cascades exhibiting a fractal-like branching structure and binary-collision-type behavior. Statistical analysis across a large number of initial directions and energies shows that no thermal spike is produced, and that the production of displacements and defects is welldescribed by the Kinchin-Pease and Norgett-Robinson-Torrens models, respectively. The simulations quantify important quantities such as the range of the primary knock-on atom and the average energy loss per collision, as well as providing a starting point for future studies of defect generation under irradiation. This information is invaluable for understanding the role of graphite under irradiation, a topic of great importance for lifetime extension of existing nuclear reactors and next-generation designs operating at high temperature.