Molecular dynamics simulation of radiation damage cascades in diamond

J. T. Buchan, M. Robinson, H. J. Christie, D. L. Roach, D. K. Ross, and N. A. Marks Department of Physics and Astronomy, Curtin University, Perth, Western Australia 6845, Australia Nanochemistry Research Institute, Curtin University, Perth, Western Australia 6845, Australia Physics and Materials Research Centre, School of Computing, Science and Engineering, University of Salford, Salford, Greater Manchester M5 4WT, United Kingdom


I. INTRODUCTION
Diamond is an important radiation-hard material with a wide range of applications, including radiation detectors for particle accelerators 1,2 and nuclear detection, 3,4 medical dosimeters, [5][6][7][8] X-ray mirrors, 9 and windows in high-power optics. 10The advantages of diamond over alternative semiconductors, in particular, silicon, stem from its unique physical and electronic properties, 11,12 including high thermal conductivity, radiation resistance, biologically compatible atomic number, large band gap, and high charge carrier mobility.These properties permit fabrication of diamond detectors with high sensitivity and low noise and capable of operating in extreme environments, such as reactors and space as well as within the human body or as a tissue equivalent; see Ref. 13 for a comprehensive overview of applications and properties.
Radiation damage in diamond has been studied for many years using ion implantation and to a lesser extent, computer simulation.One of the fundamental questions is the process by which defects are created and the vacancy concentration threshold above which diamond amorphizes and/or forms graphitic-like domains.Experimental studies of annealing in ion-implanted diamond [14][15][16][17] show that broken diamond bonds will revert back to either the sp 3 diamond or the sp 2 graphitic hybridization, and that there exists a critical dose, and hence vacancy density, N c , above which the diamond cannot be reformed.Molecular Dynamics (MD) simulations using the Tersoff potential 18 also concluded that the annealing process was dependant on N c [Ref.19].Related studies using the binary-collision SRIM package 20 similarly used vacancy densities as a starting point for a model of annealing behaviour. 21In the absence of annealing, amorphization in diamond has been shown to be driven by a similar critical threshold, 22 with experiments and MD simulation demonstrating that the diamond lattice remains largely intact down to densities as low as 2.9 g/cc, at which point $15% of the lattice sites are vacant.Applying a uniaxial strain in the simulations yielded the same threshold, demonstrating that strain is a fundamental quantity that controls amorphization.
Despite the success of computer simulation in interpreting radiation response experiments in diamond, remarkably few MD simulations have been performed on the radiation damage cascade itself.The only examples we are aware of are those of Smith 23 who modelled energetic bombardment in graphite and diamond in the 100-1000 eV energy range, and Saada et al. 19 who considered displacements arising from a 416 eV primary knock-on atom (PKA).Smith concentrated on surface ejection processes motivated by sputtering, while Saada et al. considered annealing-induced graphitisation facilitated by defects created by the PKA.While related MD studies have examined the threshold displacement energy [24][25][26] and response to swift heavy ions, 27,28 there is virtually no MD literature concerning cascade kinetics in diamond and the associated generation and annealing of defects.The importance of understanding cascade dynamics in diamond was highlighted by Prins and Derry 21 who discussed the relative importance of displacement spikes versus point defect generation.While displacement spikes are common in many metals and oxides, they are undesirable in diamond as they can be connected with residual pockets of amorphous damage which persist once the thermal spike has cooled.
In this work, we present a comprehensive MD simulation of radiation damage cascades in diamond.A large number of simulations are performed, spanning a wide number of uniformly distributed directions and energies as high as 2.5 keV.This dense data set enables statistically significant analysis of important cascade quantities, such as lifetime, spatial extent, atomic displacement, and defect recombination efficiency.We find that radiation damage cascades in diamond are strikingly different to those in most metals and 0021-8979/2015/117(24)/245901/9/$30.00 V C 2015 AIP Publishing LLC 117, 245901-1 oxides, exhibiting very short cascade lifetimes and a fractal-like trajectory structure associated with generation of isolated point defects.While there are some quantitative differences, there are many qualitative similarities with cascades in graphite, 29 a result which we found rather unexpected given the large differences in density, structure, and bonding.
Our manuscript is structured as follows.In Sec.II, we outline our systematic approach for modelling cascades and the subsequent defect analysis.Section III begins with some representative visualisations of individual cascade events, followed by a statistical treatment of kinetics and defect generation, concluding with a comparison between diamond and graphite.Section IV discusses the results in the context of previous calculations of related Group IV elements and identifies future directions for studying radiation phenomena in diamond.

II. METHODOLOGY
The calculations are performed using the environmental dependent interaction potential (EDIP) for carbon, 30,31 in combination with the Ziegler-Biersack-Littmark (ZBL) potential 32 to describe close approaches.The carbon EDIP methodology accurately models the behaviour of disordered and amorphous carbons, 33 making it well suited for describing an irradiation event.The potential is based on the earlier EDIP for silicon, 34 which contains pair and triple interactions and a spherical coordination term.The carbon variant includes aspherical coordination terms to capture dihedral rotation and long-range p-p repulsion and includes a variable short-range cut-off that ensures the sp 3 and sp 2 hybridizations reduce to the Stillinger-Weber form for integer coordination.At very small separation, the ZBL potential applies between all pairs of atoms, while for computational practicality, the three-body term in EDIP always remains active, but is comparably much smaller than the pair interaction.Due to the environment-dependence, conventional spline functions are not suitable to switch between the ZBL and EDIP pair terms, and hence Fermi-type scaling functions are used instead to smoothly interpolate and avoid points of inflexion. 29attices were equilibrated at 300 K and periodic boundary conditions were employed.Cubic supercells were always employed, with the largest cell containing 110 592 atoms and having a side length of 107.2A ˚. To confine each collision cascade within the cubic supercell, the location and initial direction of the PKA were selected such that the cascade would evolve towards the centre of the cell.A thermal layer was employed on the edge of the supercell to dissipate any excess kinetic energy.The equations of motion were integrated using the Verlet method in combination with a variable timestep algorithm developed for non-equilibrium systems. 35The metric jjF max jjDt was used to compute the timestep, which varied from 0.0005 fs during the highest energy encounters to around 0.2 fs once equilibrium was reached.Atomic motions were covered for 1 ps, which while small compared to other materials, was more than sufficient to capture all of the important dynamics in the cascade.
To acquire a sufficiently large data set for statistical sampling, a large number of cascade trajectories were sampled.PKA energies of 100, 250, 500, 750, 1000, 1500, 2000, and 2500 eV were used to initiate the radiation cascade and the initial direction of the PKA was selected from a set of points uniformly distributed on a unit sphere, formally known as the Thomson problem. 36This directional sampling technique, first used by Robinson et al. 37 to calculate threshold displacement energies in rutile TiO 2 , provides an elegant solution to the problem of selecting a set of initial velocities that are uncorrelated with the crystalline axes and representative of an arbitrary cascade.Here, we use a 25-point Thomson solution generated by our software package NanoCap 38 and shown in Fig. 1.The reduced Coulomb energy of this structure is 243.8127602,identical to that determined by Wales and Ulker 39 using a global minima search algorithm.Their data, which are available on the Cambridge Cluster Database website, 40 list the point group symmetry of the 25-point solution as C s , and hence the only symmetry operation is a mirror plane.It was particularly important to avoid high symmetry solutions, such as those with tetrahedral symmetry, as these sets of Thomson points would correlate with the lattice and significantly reduce the extent of independent sampling.
The cascade trajectories were studied using a variety of statistical techniques to quantify the response of diamond to irradiation, providing insight into energy transfer and the number of energetic atoms as a function of time.Further statistical analysis was employed to quantify the number of defects, cascade duration, and cascade length as a function of PKA energy.As in the recent work on graphite, 29 defects and displacements were determined using a vacancy radius v r of 0.9 A ˚.All quantitative data were averaged across the 25 different directions, and all time-varying data were histogrammed prior to the calculation of the mean and standard deviation.Data involving the ballistic phase only were analysed using a fixed bin width.Data which spanned the entire simulation, such as defect evolution, required two distinct bin widths, one to capture the ballistic phase and a much larger value for the annealing phase.

III. RESULTS
The results of the simulations are presented in four parts.Firstly, we show two individual cascades in diamond to provide a visual representation of a cascade.To the best of our knowledge, these are the first such calculations for diamond.In the second section, we quantify time evolution statistics by averaging over all 25 directions and determine the average cascade lifetime by analysing both the maximum kinetic energy and the number of atoms with energy exceeding fixed thresholds.The third section quantities defects and displacements as a function of time and PKA energy, while the final section quantifies the cascade size as a function of energy and compares the diamond results with those of graphite..0fs, respectively, and depict the development of the respective sub-cascades.These two frames reveal the lower sub-cascade to be more energetic with more branching evident with each successive atomic interaction.By 12.6 fs [panel (e)], the point of maximum damage is reached with a total of 20 defects.Beyond this time, the atoms do not have sufficient energy to create further defects and recombination becomes the dominant process.Recombination is extremely rapid, and after 16 fs [not shown], the number of defects has decreased from 20 to just 12. Panel (f) shows the final state of the system 1 ps after the creation of the PKA, where just 10 point defects remain.While the specifics vary according to the impact parameter, this sequence and general branching structure is typical.Intriguingly, the branching trajectories seen here are not dissimilar to those recently calculated by ourselves in graphite [see Fig. 2 in Ref. 29] The main visual contrast is that the distance between branching points is smaller in diamond compared to graphite, a difference likely due to the higher density of diamond.We will return to these similarities and differences in Sec.III D where we quantify several cascade quantities for both systems.

A. Individual cascades
A second example of a 2 keV cascade using a different initial PKA direction is shown in Fig. 3.In this simulation, the cascade nature does not contain the branching structures seen in Fig. 2, and instead involves a channelling-like process.Due to the absence of heavy collisions, the PKA travels a substantial distance, creating defects with close proximity to the central track.The upper panel of Fig. 3 shows the time of maximum defects (total of 19) which occurred 18.2 fs after the PKA was initiated.At this instant, the channellinglike behavior of the PKA has not concluded, and the cascade continues to develop afterwards.The lower panel of Fig. 3 shows the final structure once the PKA has reached the end-of-range and recombination has concluded.Comparing the two panels shows that the number of defects to have recombined is almost 50%, and as we will quantify later by averaging over many different directions, this value is typical.Comparing Fig. 3 with Fig. 2, we also see that the time of maximum defects is approximately 50% larger for the channeling case (18.2 vs 12.6 fs), indicating the intuitive result that channeling-type events spread over a longer timeframe, even though the number of final defects is similar.

B. Cascade kinetics
Having outlined some qualitative features of radiation cascades in diamond, we now turn to a more quantitative treatment.From this point onwards in the manuscript, we will discuss only statistically meaningful quantities extracted from the data for 25 different directions at each PKA energy.To gauge a sense of the strength and duration of the interactions involved in the cascade, the maximum kinetic of any atom in the system, KE max , was computed as a function of time.At each instant in time, the average value of KE max and its standard deviation was computed over all 25 trajectories.Figure 4 shows an example of this analysis for a 1 keV PKA, and also includes two specific trajectories (red and blue lines) corresponding to individual directions.These two directions were selected as they represent two extreme scenarios in the 25-trajectory data set.The blue line shows data from a localised cascade in which a heavy collision occurs as shown by the sudden drop in KE max to 270 eV at t ¼ 4.5 fs.This cascade is similar in nature to the branching calculation seen in Fig. 2. The red line shows data for a channeling-like event which is similar to that shown in Fig. 3.The signature of the channelling process in Fig. 4 is the region of small negative slope in which a fast moving atom (not necessarily the PKA) is able to move comparatively large distances without undergoing a collision.
The mean and standard deviation of KE max are shown in Fig. 4 as a solid black line and error bars, respectively.Even though the standard deviation is large, reflecting the stochastic nature of the collision process, the average value is rather smooth and exhibits a near monotonic reduction over time.
One surprising observation in this data is the rapid dissipation of kinetic energy as the cascade evolves.Only 40 fs is required for KE max to reduce to 10% of the initial PKA energy, and after 60 fs, the average value of KE max is negligible and the ballistic phase is complete.In contrast, most other materials, prototypically metals and oxides, have much longer ballistic phases, with typical durations being a substantial fraction of a picosecond and higher.It is reasonable to presume that these short lifetimes result from the dense network of strong bonds in diamond which give it a very high thermal conductivity (>2000 W m À1 K À1 ) and in a cascade context enable highly efficient energy transfer away from the PKA.
The averaged instantaneous maximum kinetic energy KE max was computed for all eight PKA energies and is shown as a function of time for three selected energies (500, 750 and 1000 eV) in Fig. 5. Also shown in the figure are black lines denoting exponential fits to the data.Each function has the form where the maximum kinetic energy at any point in time t can be described by the initial kinetic energy of the PKA and a time constant, t 0 , which quantifies the rate at which energy is dissipated.Accordingly, the averaged data for a given energy can be described by a single parameter, namely, the time constant.The goodness of fit and the simplicity of the functional form provide an elegant way to characterise energy transfer in the cascade and may well be useful in radiation damage studies on other materials.Data for the other five energies were also accurately described by this functional form, and via the fitting procedure, a value of t 0 was extracted for each PKA energy.The values of t 0 are shown as red circles in Fig. 6 where they are plotted as a function of energy.Also shown in the figure is a power-law function as previously used to analyze cascades in graphite.The power-law has the form aE x , and yields an exponent x for 0.36 when fitted to the time constant values.Although there is some scatter at the higher energies, the fit at low energies is excellent.We suspect that this scatter reflects sampling noise, since at the higher energies some of the cascades (1 at 1 keV, 2 at 1.5 keV, 4 at 2 keV and 1 at 2.5 keV) contacted the cell boundary and were removed from analysis; for all other energies, all 25 directions remained fully contained.Removing these cascades tends to reduce the time constant, since the trajectories which are omitted have large spatial extents and lose energy more slowly (see Figs. 3 and 5).Significantly, the value of the exponent is much smaller than the value of 2 3 corresponding to cooling of a thermal spike; see Ref. 41 for a derivation based on the heat-diffusion equation.The absence of a thermal spike (or equivalently, cascade evolution via pointdefect creation), is again likely to be related to the high thermal conductivity of diamond since the local temperature never rises high enough to create a molten cascade core.We will return to this property in Sec.III D where we compare the cascade lifetime in diamond with graphite using a slightly different metric.
To quantify the rapid time-evolution of the radiation cascades, the number of atoms exceeding pre-defined threshold values of 1 and 10 eV was counted and averaged across all directions.The higher threshold provides a measure of fast moving atoms associated with defect production and atoms moving between lattice sites, whereas the lower threshold conveys a sense of heat dissipation to the surrounding atoms in the lattice.These categories can be considered a measure of "fast atoms" and "warm atoms," respectively.Data for the 1 keV set of directions are shown in Fig. 7, with fast atoms shown in red (lower curve) and warm atoms in blue.As before, error bars denote one standard deviation.The short duration of fast atoms mirrors the rapid exponential decay seen in Fig. 5, and the number of these atoms is rather small, reaching a maximum value of six.In contrast, the number of warm atoms becomes relatively large, around 100, before gradually reducing to zero around 0.25 ps.A high degree of statistical variability is apparent towards the end of the ballistic phase when the number of warm atoms is at its maximum (the 1r confidence level extends from around 75 to 130), but from around 0.05 ps onwards, the standard deviation is much smaller, and the dominant process is thermal diffusion in which the residual kinetic energy is conducted away into the surrounding matrix.
In the earlier graphite study, we extracted a useful quantity t max from individual trajectories by searching for the time at which the number of fast atoms was maximal; this time quantifies growth of the cascade during the ballistic phase.One drawback of the previous approach was that it suffered from large statistical noise as the number of fast atoms was small and assumes an integer value for single trajectories.With the present data, however, it is possible to define t max much more precisely, since the number of fast atoms is now comparatively smooth and can take fractional values.An example of this analysis is shown in Fig. 8 where data are shown for three different PKA energies.Also shown in the figure are parabolically-shaped dotted lines determined by fitting functions of the form (2) where a, t max , and N max are fitting parameters.Using this approach, t max can be easily extracted as a function of PKA and interpreted to understand the behaviour of the ballistic phase of the cascade.We will return to this discussion in Sec.III D when we compare behavior in diamond with that in graphite.

C. Defect analysis
Structural analysis of the instantaneous positions of the atoms provides a complementary view to the metrics based FIG. 6.Time constant, t 0 , as a function of PKA energy determined by fitting Eq. ( 1) to maximum kinetic energy data, such as shown in Fig. 5.The solid line shows a power-law fit as described in the text.on kinetic energy.Two of the most important quantities are the number of defects, defined as atoms more than v r ¼ 0.9 A from a lattice site, and the number of displacements, defined as atoms which have moved more than v r from their original position.During a cascade, atoms can migrate from one lattice site to another, and hence the number of defects is a lower bound to the number of displacements. Figure 9 shows the time evolution of displacements and defects for cascades with an initial PKA energy of 2 keV.Once again, data are averaged across many directions (here 21, since four cascades interacted with the walls), and error bars denote one standard deviation.Note that the finer spacing between the data points at small times reflects the smaller bin width required to average across the variable integration time step and accurately capture the ballistic phase of the cascade.Early in the cascade, the average number of defects and displacements are equal, but by the time the maximum number of defects is attained at 0.05 ps, the number of displacements is already significantly larger, by $15%.Beyond this point, the number of displacements is close to constant, except for a small reduction in displacements in which mobile atoms return to their original lattice.Due to recombination of vacancies and interstitials, the number of defects reduces substantially as the cascade evolves, eventually falling to just under half the maximum value.Although the averages are welldefined and smooth functions in time, the magnitude of the standard deviation is substantial.Without the averaging procedure employed here, it would be impossible to extract quantitatively useful information such as seen in Fig. 9.
The analysis shown in Fig. 9 was performed for all of the PKA energies and the number of displacements and defects at the conclusion of the cascade was determined.These data are plotted in Fig. 10, where it is apparent that the number of displacements and defects each exhibit a strong linear dependence on the PKA energy; the black lines are fits to the data, while the error bars are one standard deviation.The linear trend is extremely strong, with the lines passing through all of the data points; indeed, if the standard error in the mean is plotted, the error bars become so small as to be barely distinguishable from the data point itself.These two data sets are reminiscent of the historical Kinchin-Pease 42 (KP) and Norgett-Robinson-Thomson 43 (NRT) models in which the threshold displacement energy (E d ) is the key quantity relating the PKA energy to the number of displacements and defects, respectively.Assuming a KP functional form yields a value of E d of 41 eV, while the NRT functional form yields a rather larger value of 69 eV.While neither of these values should be considered definitive (since explicit trajectories have not been used to calculated defect formation probabilities as in Ref. 37), they are consistent with experimental and computational values in the literature [24][25][26][44][45][46][47] which span a broad range extending from 30 to 80 eV.

D. Comparison to graphite
To place the present results in a broader context, it is instructive to compare some of the properties of the diamond cascades with related simulations recently performed by ourselves on graphite. 29We first consider the cascade length and PKA range as shown in Fig. 11.The size of the cascade is defined as the largest distance between any two defects in the cascade, while the PKA range is determined by calculating the difference between the initial and final positions of the PKA.Each data point in Fig. 11 represents an average across 25 directions and the error bars indicate the standard deviation.The PKA range and the cascade length are each well-described by straight-line fits passing through the origin, represented by the black lines.For the PKA range, the gradient of the linear fit for diamond is 18 A ˚/keV, approximately 40% smaller than the corresponding value of 30 A ˚/keV for graphite. 48Similar behaviour occurs for the cascade length, where the gradient for diamond is 25 A ˚/keV, around one-third smaller than the graphite value of 38 A ˚/keV.As noted earlier, the quantity t max provides a useful measure of the timescale of the ballistic phase.Figure 12 compares a power-law fit extracted from our previous graphite study with corresponding data for diamond extracted using Eq. ( 2) and illustrated in Fig. 8.The data for diamond are also well-described by a power-law expression, and the exponents x are very similar, 0.41 for diamond versus 0.37 for graphite.Since the values of x are so close, the ratio between the diamond and graphite values does not vary much with PKA energy.Examination of the two data sets shows that t max in diamond is approximately 45% shorter than in graphite, consistent with the more compact cascade dimensions seen in Fig. 11.
While we do not have a complete explanation for the power-law behaviour, one plausible hypothesis is that it arises from the branching nature of the cascades.Applications as diverse as internet connectivity, 49 ecology, 50 and metallic glasses 51 show that fractal behaviour is associated with powerlaw relationships.It therefore seems reasonable to propose that the power-law behaviour evident in the cascade timescale data is associated with the fractal-like behaviour of the cascade trajectories themselves.Given this potential mathematical connection, the behaviour seen for graphite and diamond is worthy of further investigation, particularly in other similar materials to test if there is some underlying universality.
Figure 13 plots the ratio of defects to displacements as a function of PKA energy for diamond and graphite.In the original NRT article, 43 this ratio is referred to as the "displacement efficiency," but we prefer the label "residual defect factor," since this emphasises that a small ratio implies few residual defects, and vice-versa.According to the NRT model, this factor is independent of PKA energy, material and temperature, and has the value 0.8, corresponding to 20% recombination.Figure 13 shows that the first of these assumptions (no PKA energy dependence) holds true for diamond and graphite, but the second assumption is clearly incorrect.Across a wide energy range, we find the residual defect factor is 0.74 for graphite, not dissimilar to the NRT value, and just 0.50 for diamond.This implies around one-quarter of the displaced atoms in graphite recover to find new lattice sites, while recombination in diamond is much more significant, with half of the displaced atoms finding new lattice sites.The strong propensity in diamond for recrystallization is further magnified by a comparatively low rate for generating atomic displacements in the first place.As seen earlier, using the Kinchin-Pease model (Eq.( 3)) as a simple metric to quantify, the displacement efficiency yields a value of E d of 41 eV for diamond, considerably higher than the equivalent value for graphite which is around 25 eV (see Fig. 9 in Ref. 29).Accordingly, around 40% fewer displacements are generated in diamond relative to graphite for a given PKA energy, and together with the higher 12.Time of maximum number of fast atoms, t max , as a function of PKA energy.The black line is a power-law fit to the diamond data, shown as red circles.The blue line shows the power-law relationship for the same quantity from our recent graphite studies. 29G.13.The residual defect factor, defined as the ratio of defects to displacements, as a function of PKA energy.Data for graphite and diamond are shown as blue squares and red circles, respectively.The graphite data set is taken from our previous MD simulations 29 as discussed in the text.recombination efficiency, diamond overall yields 60% fewer defects than graphite.

IV. DISCUSSION AND CONCLUSION
The predominant insight to emerge from this study is that radiation cascades in diamond are distinctive as compared to many other materials, in the sense that the properties of the radiation damage process fall at the far end of the spectrum.While there are clear quantitative differences between damage in diamond and graphite, in many ways these two principal carbon allotropes exhibit similar radiation response and can be grouped in a relatively small class of materials in which the cascade response involves point defects and isolated subcascades.The most important factors common to diamond and graphite cascades are the branching structures of the trajectories, isolated point defect production, an extremely short ballistic phase and the absence of melting as per a thermal spike.All four of these attributes are fundamentally connected, arising from efficient transfer of kinetic energy away from the collision site.As a result, the binary collision approximation (BCA) is a useful starting point for studies at higher energies, a result we previously discussed at some length in our study of graphite cascades. 29While the BCA is commonly (and correctly) assumed in discussion of graphite, in particular nuclear graphite, we are not aware of a similar understanding in the literature for diamond.This finding is an important theoretical result for a diamond, where atomistic knowledge of radiation effects is surprisingly modest, with the vast majority of studies concerned with either fabrication or characterisation in the laboratory or high-energy simulation, such as with GEANT4. 52As for why diamond and graphite exhibit this BCA-type behavior, a full explanation is still lacking, but molecular dynamics is the ideal tool to elucidate the specific characteristics responsible by exploring other carbon phases, such as nanotubes, glassy carbon, and amorphous carbon.
Comparison of the present diamond results with earlier cascade studies [53][54][55] in silicon and germanium is instructive, as the latter two materials have the same diamond crystal structure but weaker covalent bonds.In both silicon and germanium, the cascade evolves in a manner similar to many materials, producing an extended volume of disorder that results in an amorphous pocket when the cascade cools.This observation that isolated point defect production is not common in silicon and germanium shows that it is strong carbon-carbon bonds in diamond which are responsible for the BCA-like behaviour, rather than the crystal structure itself.Note that this cascade evolution in diamond is fundamentally different to radiation resistant oxides, such as spinels and pyrochlores where the origin of similar behaviour derives from structural factors, such as cation anti-sites and structural vacancies.
Despite the many conceptual similarities between diamond and graphite, the quantitative details are quite different.As we have discussed separately, 29 compared to most materials, graphite can be considered rather unusual due its very fast cascade lifetime, fractal-like trajectory structure, and absence of a thermal spike.Furthermore, graphite resists radiation damage extremely well, which is why it was the original nuclear material, being employed as a moderator, reflector, and structural material.Against this backdrop, diamond is even more radiation resistant than graphite, exhibiting ballistic phases which are 45% faster (Fig. 12), cascade dimensions that are $35% smaller (Fig. 11), and defect production which is 60% lower due to more efficient recombination (Fig. 13) and a higher Kinchin-Pease threshold displacement energy (Fig. 10).All of these differences combine to the high radiation hardness of diamond and contribute to the widespread interest in exploiting diamond as a detector, dosimeter, and other radiation-hard situations.
Looking beyond the particular case of diamond, one of the technical achievements in this work is the use of directional sampling and averaging to extract robust statistical trends.The averaging of the maximum kinetic energy as a function of time is a case in point: even though an individual cascade is highly specific to the direction of the PKA (Fig. 4), averaging over a sufficient number of trajectories reveals a strong exponential trend (Fig. 5).This trend could not have been anticipated by inspection of individual cascades due to the large fluctuations associated with heavy collisions.It would be intriguing to learn whether other materials important in a radiation damage context also exhibit this exponential energy dissipation process.If that proved to be the case, this averaging approach would provide a simple metric to compare different materials systems.A similar logic applies to the other averaged quantities presented above.The number of fast and warm atoms provides insight into the cascade that is intuitive to interpret, as is the time when the number of fast atoms is highest.Although not presented here, the time at which the number of defects is maximal is another useful quantity that can be easily extracted since the averaged data are so smooth (Fig. 9).Many other useful quantities undoubtedly be defined, and providing this quantitation in a comprehensive way across multiple materials systems would be invaluable for understanding radiation damage in solids in general, and diamond in particular.
In summary, we have performed a comprehensive molecular dynamics study of radiation damage cascades in diamond.Similar to graphite, cascades in diamond generate isolated point defects and a branched trajectory structure created by distinct heavy collisions.Statistical averaging over a uniformly distributed set of initial directions was applied to quantify a variety of diamond cascade characteristics, including energy dissipation (decays exponentially with time), cascade lifetime (power-law dependence on PKA energy), spatial extent (linear dependence on PKA energy), and defect formation/recombination (also linearly dependent on energy).Quantitative comparisons with cascades in graphite show that diamond is an even more extreme material, with briefer cascades and fewer defects.This high degree of radiation hardness is one of the key reasons diamond is such an attractive material for dosimetry and detectors, and motivates further computational studies to pinpoint the threshold displacement energy and develop a complete atomistic picture of radiation damage in diamond.
Australian Government and the Government of Western Australia.NAM also thanks the Australian Research Council for support (FT120100924).

FIG. 1 .
FIG. 1. Solution to the Thomson Problem for N ¼ 25 in which point charges (shown in red) are uniformly distributed on the unit sphere.In this work, the initial direction of the PKA corresponds to one of these 25 points.

Figure 2
Figure 2 shows an example of a cascade in diamond produced by a 2 keV PKA.Interstitials and vacancies are represented by blue circles and red squares, respectively, and the black lines show the trajectory of the defect atoms, providing a sense of the evolution and spatial extent of the cascade.All other atoms which remain on their lattice site are not shown.Panel (a) shows the cascade 1.8 fs after the initial creation of the PKA, at which point the PKA has moved a small distance

FIG. 2 .
FIG. 2. Time evolution of a typical cascade in diamond for a 2 keV PKA.Interstitials are denoted with blue circles and vacancies shown as red squares.Atoms that remain on their lattice sites are not shown.Panel (a) depicts the initial stages of the PKA event and the final frame in panel (f) indicates the stable defects that remain once recombination has occurred.

FIG. 4 .
FIG.4.Maximum kinetic energy as function of time for cascades with an initial PKA energy of 1 keV.Blue and red lines present data for two specific cascades as described in the text.The solid black line and error bars denote the mean and standard deviation, respectively, averaged over 25 uniformly distributed initial directions.

FIG. 7 .
FIG. 7. Number of atoms with kinetic energies greater than 1 and 10 eV as a function of time for a PKA energy of 1 keV.Data averaged across 25 directions with error bars indicating one standard deviation.

FIG. 9 .
FIG. 9. Number of displacements (red; upper curve) and defects (blue; lower curve) as a function of time for a PKA energy of 2 keV.Data are averaged across 21 trajectories using two distinct bin widths.Error bars denote one standard deviation.

FIG. 11 .
FIG. 11. (a)The PKA range and (b) cascade length as a function of PKA energy.The error bars denote one standard deviation.The linear fit for diamond (black) is compared to another fit for graphite (blue) from previous MD simulation studies.29