Molecular Dynamics Analysis of Collison Cascade in Graphite: Insights from Multiple Irradiation Scenarios at Low Temperature

Abstract

In our study, we utilize molecular dynamics simulations, specifically through the Reactive Empirical Bond Order, to unravel atomic-scale dynamics in graphite, an essential component in many advanced technologies, under varying irradiation scenarios. We shed light on the behavior of graphite when exposed to Primary Knock-on Atom (PKA) energies of 10, 20, 40, and 80 keV. The findings highlight the radiation vulnerability of graphite, especially when influenced by hydride inclusion. Both pristine graphite and its hydride variant exhibited an increase in Frenkel pairs (FPs) with escalating PKA energies. Notably, carbon PKAs manifested significant FP effects, whereas hydrogen PKAs influenced defect formation through variable diffusivity. In tested radiation scenarios, particularly in Mode C and the R1 region, cascade patterns identified distinct defect forms of diamond-like and elongated-diamond-like shapes, distinct from the typical PKA collision clusters. Furthermore, our cascade findings emphasize the formation of three-coordinated graphite rings, particularly as PKA energies increase. The graphite population statistics reveal a decline in threefold-coordinated atoms and an increase in other types of defects, with 7-carbon atom rings being the most common. Our research highlights the significance of understanding three-coordinated graphite rings, especially as PKA energies rise. Graphite population statistics reveal a decline in threefold-coordinated atoms and a rise in other defects. Notably, 7-carbon atom rings are the most common. From a clustering perspective, self-interstitial atom (SIA) clusters predominated in pristine graphite, while this trend balanced in the hydride variant. Our research highlights the importance of understanding atomic behaviors in graphite under several radiation scenarios. This knowledge is needed for advancing reliable and efficient technological applications, particularly in the field of nuclear technology. Our research underscores the need to understand atomic behaviors in graphite under radiation, paving the way for detailed study on reliable efficient technological applications.

1. Introduction

Graphite has always been at the forefront of materials suitable for nuclear applications, owing to its unique physicochemical properties. Its legacy dates to 1942 with the inception of the Chicago Pile 1, and it has been considered for advanced Generation IV reactor concepts [1]. Such persistent reliance on graphite is attributed to its remarkable radiation responses, which induce significant changes in its thermal, mechanical, and creep properties [2]. Neutron-induced interstitials in graphite can raise its energy levels, leading to phenomena like Wigner energy release if not effectively managed.

Graphite’s behavior under irradiation, especially the intriguing radiation-induced dimensional change involving shrinkage, swelling, and cross-over at elevated damage levels, is complex yet crucial for reactor performance [3]. Pre-existing Mrozowski cracks and porosity in reactor-grade graphite interact with defects, causing atomic-scale changes in structure [4].

Traditional methods, like Scanning Tunneling Microscopy and Rutherford backscattering, have been employed for investigating atomic-scale radiation damage in carbon [5]. STM, although adept at surface characterization, offers limited scope in area analysis. Similarly, Rutherford backscattering often struggles with bulk materials due to the inherent challenges posed by polycrystalline structures like HOPG. In this landscape, Raman spectroscopy has emerged as a promising technique. It indirectly pinpoints defects in carbon materials, with the D to G Raman mode intensities being indicative of defect populations. The modified Tuinstra–Koenig (TK) model has further elucidated this, revealing that at higher damage levels, the D mode intensity decreases as defects become more proximate to the sixfold ring size.

Yet, it is the advent of Molecular Dynamics (MD) simulations that has truly revolutionized our understanding; offering unparalleled insights at the atomic level, MD simulations highlight pivotal factors like primary knock-on atom (PKA) energy [6,7,8,9,10]. Nevertheless, the selection of an appropriate interatomic potential remains contentious [11]. Historical attempts, such as using the Tersoff potential, have provided valuable findings on the sputtering of graphite [12]. Later efforts incorporated potentials like the Brenner, renowned for more properly capturing the C60 structure [13]. However, these potentials often falter when faced with long-range interactions. A significant breakthrough was the modified Tersoff potential with a long-range extension, shedding light on previously unnoticed hillock formations on graphite surfaces [14].

Despite these advancements, challenges persist. For instance, while density function theory (DFT) has offered a glimpse into early-stage defect evolution, it grapples with the intricacies of weak van der Waals interactions [15]. Yet, with the introduction of the AIREBO potential, accounting for both inter-layer covalent bonding and van der Waals forces, there is renewed hope [16]. However, values like threshold displacement energy were successfully evaluated with MD, with values spanning from 10 eV to 70 eV, as discussed by Banhart and Zinkle [17] As we move forward, understanding irradiation schemes on graphite, especially in high-energy cascades, becomes paramount.

While graphite, in its pure form (i.e., pristine), has been the subject of extensive studies, irradiation-induced impurity effects such as the hydrogen production in reactors can lead to its accumulation within graphite, morphing it into a graphite hydride [18]. Hydrogen can influence graphite’s thermal stability, mechanical resilience, and radiation resistance [19]. With the critical roles of graphite in reactors like gas-cooled reactors and very high-temperature reactors (VHTRs), any change in its structure and properties, due to impurities or otherwise, is worthy of consideration. Furthermore, variant effects such as irradiation patterns can be still searchable with MD simulations.

In this work, we extend the ongoing discourse by leveraging MD simulations to closely examine radiation-induced collision cascades in graphite. Our study spans a range of primary knock-on atom (PKA) energies, by varying displacement directions, and even delves into multi-PKA irradiation scenarios for a more realistic representation.

In our study, we conducted a series of MD simulations to explore the collision cascades in both pristine graphite and graphite hydride. The investigation was carried out over a range of PKA energies, specifically at 10 keV, 20 keV, 40 keV, and 80 keV, to provide a comprehensive understanding of the energy-dependent response of the materials. To further deepen our insights, we not only examined the effect of the collision displacement direction and vibration effects but also considered different scenarios of irradiation. In addition to the traditional single PKA approach, we explored complex irradiation scenarios where the same total amount of PKA energy was distributed among several PKAs at different locations within the graphite structure, simultaneously. This multi-PKA approach allowed us to mimic more realistic irradiation conditions, potentially offering a more nuanced view of how graphite responds to radiation. By systematically varying these parameters, our work paints a wider picture of the irradiation behavior of these materials, contributing valuable insights that are relevant to their application in nuclear environments.

2. Simulation Methods

In this research, molecular dynamics (MD) simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package [20].

The carbon interatomic potential was described by the modified Adaptive Intermolecular Reactive Empirical Bond Order potential (AIREBO) and used to characterize the interatomic interactions for graphite and graphite hydride [21], which was developed from the original AIREBO potential by J. Stuart [16]. The energy potential parameters for carbon–carbon (C-C) formation and migration, as well as hydrogen–hydrogen (H-H) interactions, were established through alignment with a mix of experimental data and first-principles calculations. Additionally, the cross-interaction parameters for carbon–hydrogen (C-H) were refined to conform to a specific set of first-principles data.

We conducted MD simulations by subjecting a simulation cell structure to initial conditions of 30 K and 0 Pa in order to equilibrate it. For the recoil simulations, the starting temperature was established at 30 K, but for some simulations, it was reduced to 0 K to obtain observations unaffected by the thermal movements of atoms. These simulations were categorized as the 30 K recoil simulations and the 0 K simulations, respectively.

The objective of conducting simulations at 30 K, termed 30 K recoil simulations, was to align the simulations closely with experimental values noted at low temperatures [22]. To mitigate the impact of atomic vibrations on the accuracy of our collision cascade predictions at this temperature [23], we performed several iterations for each directional displacement. This entailed conducting four distinct simulations for assessing the collision cascade in every direction, with the time gap between successive recoil incidents fixed at 50 fs. These simulations commenced at varying instances: 0, 50, 100, and 150 fs [23]. Averaging the outcomes from these runs helped in reducing the effects of thermal vibrations and temporal correlations on the cascade data.

The pure graphite structure, formed by replication rectangular graphite unit cell shown in (Figure S1a), is constructed by stacking graphene layers placed in periodic boundary conditions at all boundaries in a cell size measuring 116 Å × 118 Å. Approximately 40 of these layers are stacked together, resulting in a final structure dimension of 116 Å × 200 Å × 118 Å. This assembly comprises approximately 400,000 carbon atoms. In this configuration, the hydrogen atom substitutes for a carbon atom in each graphene unit cell, leading to the formation of C-H of the hydride structure.

The H atoms within the graphite structure were stabilized and optimized using the DFT calculation with Quantum Espresso-implemented Materials-Square Platform [24]. The selected unit cell after DFT calculation was added to the Supplementary Materials as the most stable structure (Figure S1b). The PKA was predominantly selected to be a carbon atom, with hydrogen atoms used in specific instances; both were positioned at the simulation box’s center. To evaluate the collision cascade, we averaged the defect evolution measurements across the chosen dataset. The impact of vibrations was considered by testing four different vibration timings for each direction of displacement, ensuring the statistical reliability of our findings [23]. The investigation included conducting simulations for both pristine graphite and graphite hydride systems, with results compiled from 20 simulations for each structure. The simulation cell size was determined to maintain the system’s average temperature below 200–300 K post-collision, preventing the displacement cascades from extending beyond the cell boundaries, even at the highest PKA energies.

An adaptive timestep was employed, setting the upper limit for displacement per step (xmax) at 0.01 Å and the maximum timestep (tmax) at 0.02 ps. These parameters proved to be effective for evaluating defect formation across various xmax and tmax values [25]. At the start of each MD simulation of recoil events, the PKA was imparted with recoil energy through its velocity components. Defect formation was assessed using Voronoi analysis [26]: a site was deemed to contain a self-interstitial atom (SIA) if it held more than one atom, considered a vacancy if it was empty, and regarded as undamaged if neither condition was met, for each level of PKA energy. Throughout the displacement cascade triggered by irradiation, the system’s temperature, volume, and number of atoms were kept constant (NVT ensemble) over 30 ps to ensure the observation of all three ballistic phases was feasible. To calculate the number of atomic displacements leading to defect creation in materials under irradiation, the Norgett–Robinson–Torrens (NRT) formula is widely utilized [27].

$$\mathrm{NRT}\mathrm{displacement}=\frac{0.8E_{pka}}{2E_d}$$

(1)

Here, E_pka represents the energy deposited in the material due to nuclear interactions during a collision, which is approximately equal to the energy possessed by PKA as it initiates a cascade of secondary atom collisions. E_d signifies the material’s threshold displacement energy (TDE), identified as 35 eV for pure pristine graphite, derived from applying the AIREBO potential within the computational framework outlined in the referenced study [28].

Table 1. Outline of the parameters used to simulate collision cascades. A structure measuring 116 Å by 200 Å by 118 Å was constructed. To cover various scenarios, sixteen unique simulation scenarios were devised, encompassing four vibration timing options and four displacement directions. The energy levels examined ranged between 10 keV and 80 keV.

PKA Energy (keV)	Number of Atoms in Structure	Number of Independent MD Runs	Displacements Counts (NRT)	Simulation Duration (ps)
10	400,000	20	114	30
20	400,000	20	228	30
40	400,000	20	457	30
80	400,000	20	914	30

details the modeling parameters utilized within this research, including a summary of the standard NRT displacement calculations. However, the NRT approach presumes a uniform impact of atomic displacements on FP’s formation, an assumption that oversimplifies the actual process. The formation of FPs is influenced by numerous variables, including the displacement’s direction and energy, the types of atoms involved, and the material’s local structure. Consequently, the NRT model might not accurately reflect the quantity of FPs generated, as it overlooks these critical factors. In contrast, MD simulations offer an intricate and accurate method for examining FP formation and dynamics under various irradiation scenarios, aiming to accurately quantify FP formation. The analysis of the final atomic arrangements was conducted using OVITO, a tool designed for visualization and post-processing [29]. The characteristics and dimensions of the cascade were identified using the analytical features available in OVITO.

3. Results and Discussion

3.1. Collison Cascade Dynamics

3.1.1. Frenkel Pairs (FPs) Evolution under Different PKA Energies for Graphite

This study utilized the standard error of the mean (SEM) to quantify the variations in the creation of Frenkel Pairs (FPs). The SEM was computed as the standard deviation of the average number of FPs from 20 simulations, divided by the square root of n. The damage evolution patterns (Figure 1) for both (a) pristine graphite and (b) graphite hydride at PKA energies of 10, 20, 40, and 80 keV. The patterns are quite consistent across these energy levels; higher PKA energies yielded higher peaks and prolonged times to stable recombination, whereas lower PKA energies led to shorter peaks and faster recombination. Approximately 60% of the maximum FPs survived in the peak across the range of the evaluated PKA energies. During the ballistic phase, the kinetic energy of the PKA triggers a rapid rise in atomic displacements, reaching the highest concentration of defects before most revert to their original lattice positions in a process known as recombination or stabilization. Interestingly, the peak is narrower in pristine and wider in the hydride structures. This behavior can be attributed to the influence of hydrogen in the graphite lattice. In pristine, the atomic interactions and bonds are more consistent, leading to a quicker recombination resulting in a narrower peak. Hydrogen atoms can affect the carbon–carbon bonds in the graphite lattice, which are attributed to hydrogen’s ability to facilitate bond breaking, potentially hindering the recombination process and leading to a wider peak. Also, the lighter mass of H compared to carbon atoms leads to more structural deformation.

Different PKA energy levels result in varying cascade core volumes and defect amount counts. Higher PKA energies generate results in more primary collision events and increase the probability of secondary PKAs and subsequent mini cascades [30]. This also can contribute a contribution to the broader peak observed in hydride compared to the narrower peak in pristine. It explains why an 80 keV PKA resulted in amplified structural failure in graphite hydride but remained stable and had lower counts in pristine graphite.

In Figure 1b, the occurrence of two distinct peaks (at 7 ps and 28 ps) at 20 keV requires a closer examination of the defect and structural phases. These peaks arose from interactions mainly with H atoms, leading to more FP. Though these were most common in only one displacement direction, the presented results are averages. Surprisingly, despite fewer H than C atoms, this behavior underscores the unique aspects of the radiation damage environment. In hydride, the FP counts were 6–7 times higher than in pristine, attributed to hydrogen’s role in reducing the TDE to below 30 eV, as has been concluded in other studies [31]. This increased defect rate suggests the hydride’s reduced stability compared to pristine. Additionally, hydrogen atoms can penetrate the graphene structure when exceeding 200 eV [32], whereas other recent studies have stated that the penetration threshold energy can be penetrated at much lower energies [33]. This information is important in irradiated graphene where passivating the carbon dangling bonds with hydrogen atoms in irradiated graphene delocalizes the charge density, reduces the density of states at the Fermi level, and increases the band dispersion, thereby tuning the properties of graphene for additional electronic functionality [34].

3.1.2. The Impact of Displacement Direction on the Frenkel Pairs Kinetics

Primary radiation damage is influenced by the crystal structure and unpredictable collision events that form sub-cascades and channels [35]. The stochastic nature of the cascade mechanism makes predicting PKA energy distribution tough, so more statistics are needed. To study the effects of PKA energy and direction on radiation damage, 20 simulations were run for each PKA energy. We looked at the number of remaining FPs or defects that lasted past the relaxation time, as shown in Figure 2.

We looked at displacement directions [1000], [1100], [1120], [1123], and [0001]. Usually, the defect count was found to follow the order: [0001] > [1000] > [1100] > [1121] > [1123]. The direction [1123] had the lowest defect count at energies like 10 keV, 20 keV, and 40 keV. At an energy level of 80 keV, the [1121] orientation exhibited a higher number of defects, potentially as a result of the “concentrated collision sequence” observed at elevated energies. In such a sequence, atomic collisions predominantly occur along a single direction, with the sequence’s pattern showing minimal sensitivity to changes in atomic density. This feature elucidates the observed variation in Frenkel pairs (FPs) in the pristine [1121] orientation in relation to the energy of the primary knock-on atom (PKA).

Head-on collisions, like in the [0001] direction, often result in greater energy deposition. This direction typically has the lowest TDE energy in many crystals, including graphite and graphene. Consequently, a lower TDE means atoms are more readily displaced, as indicated in references [36,37]. As PKA energy rises, differences in defect formation between crystal directions become clearer, as seen in Figure 2a,b. For example, at 80 keV, there is a notable gap in defect counts among the five displacement directions, but this gap narrows at 10 keV. This is not unique to carbon; materials like Iron and Tungsten show similar patterns [38,39,40]. The amplified differentiation at higher energies might arise from inherent directional threshold variations and the effects of focused cascade collisions. In hydride, radiation causes many more defects than in pristine. Since hydrogen atoms are smaller than carbon atoms, they can fit between carbon spaces, disrupting the structure and creating more defects due to the weakened carbon–carbon bonds [41]. Our results also observed that, under the same PKA energy and collision cascade conditions, the number of broken C-H bonds was higher than that of C-C bonds (Video S0a,b). In hydride, the effect of direction on defect formation is clearer. For example, in the <1123> direction (i.e., d direction), the material remains more stable. Hydrogen atoms mainly stay within their layer, reducing the risk of hydrogen buildup. During irradiation, these atoms tend to move away from vacancies, leading to hydrogen voids. They also cluster, which can weaken the material. The presence of hydrogen, due to its weaker bond with carbon, further enhances the chances of void creation. At 40 keV PKA energy, the [1100] direction behaves differently. Irradiation here results in 700 FPs, much more than about 250 FPs in other directions (Figure 2II). This is due to the densely packed atoms in the [1100] direction, which is also graphite’s c-axis where layers of graphene stack. The c-axis’s weaker bonds make it more prone to radiation damage. Graphite’s structure, with loosely connected graphene layers, shows layer-by-layer damage upon radiation. The top layers take most of the hit, while deeper layers, with diminishing radiation energy, have fewer defects. These layers can also slide, causing radiation to affect neighboring layers [42]. However, at 80 keV PKA energy, radiation damages multiple graphene layers, resulting in an average of 1200 to 1700 FPs across all PKA directions (Figure 2II). But, thanks to the natural strength of graphene and how radiation spreads, deeper layers might mostly stay untouched if the energy fades before reaching them.

3.1.3. Effects of Hydrogen Atoms as PKAs on Cascade Dynamics within the Hydride Structure

It is noteworthy that the employment of H atoms as PKA yields significantly fewer defects compared to C atoms. In some cases of displacement in the [1123] direction, almost no defects are observed, as illustrated in Figure S2. The trajectory of H atoms, as they navigate through the layers and rings of the CH hydride structure, varied based on the displacement direction. Complex directions presented short, twisted paths, while the more direct trajectories like the [0001] direction were linear and extended. This observation is pivotal because it accentuates the idea that the identity of the colliding atom—be it H or C—can profoundly influence the cascade and resultant damage in the hydride structure. Validating that the chances of the H or C atom being a PKA is a random and instant action in the chaotic irradiation environment. Thus, our judgment on the damage observed in Figure 1b should not be always the case.

Different displacement directions uniquely influence the PKA pathways in hydride. For example, the [1123] direction predominantly presents shorter distances between exit and re-entry points. However, the [1120] and [1000] directions typically show deeper penetration into the graphite layers. This highlights the PKA direction’s impact on radiation damage severity, in accordance with studies highlighting the influence of PKA type and direction on defect formation in irradiated materials [43,44,45]. H atoms, being lighter than carbon in the hydride structure, inherently have lower momentum [46,47,48]. This makes them less likely to displace carbon atoms within the lattice. Their superior mobility, compared to carbon atoms, enhances their diffusion ability within the structure, promoting recombination with vacancies [49]. Additionally, H atoms can introduce SIA defects, which, compared to FP, exhibit less stability and a greater recombination tendency.

Notably, as we have tested, large energies of PKA (H atom) moving from 10 keV to 80 keV do not always correlate with increased FPs. Instead, the predominant factor influencing diffusion patterns mainly appears to be the displacement direction. To clarify, we present the [0001] direction as a representative case (Figure S2). Although this direction displays a relatively linear trajectory, other directions exhibit complex pathways. The visual representation, both in wrapped and unwrapped trajectory formats, highlights the defect outcomes. Videos S1–S5 further elucidate these trajectories. The H trajectory elucidates the three phases of hydrogen interaction (Figure S2d), where the initial energetic movement transitions to a zigzag pattern before the atom becomes trapped between graphene layers. In this trapped state, the hydrogen atom can marginally alter the C-C bond length before the system reverts to its original state (Figure S3e).

The timeframes for H trajectories travels across various directions have been observed as 27.3 ps ([1000]), 33.3 ps ([1100]), 24 ps ([1120]), 24 ps ([1123]), and 40 ps ([0001]). These timings indicate that more intricate directions are prone to frequent collisions, resulting in short trajectories. Conversely, the [0001] direction, which showed uninterrupted movement, hints at enhanced internal movement potential within the structure’s layers.

3.2. Non-Conventional Collision Cascade Scenarios

To gain a comprehensive understanding of radiation damage in graphite structures, it is pivotal to model the unpredictable nature of irradiation environments more closely. Heggie et al. [50] have pointed out that graphite’s layered structure offers unique radiation interaction dynamics when exposed to high doses of energetic particles. In our study, to delve deeper, we explored a variety of collision cascade scenarios using a PKA energy of 80 keV in the [0001] direction. This approach is a departure from traditional methods that trigger the cascade event with just a single PKA. Our method contrasts this by referencing conventional FP outcomes to state differences in radiation damage. We considered four distinct energy deposition modes in the graphite structure:

• Mode D: A conventional single 80 keV PKA, representing the well-understood collision pattern.
• Mode C: Dual PKAs, each with 40 keV.
• Mode B: Quadruple PKAs, each carrying 20 keV.
• Mode A: Eight PKAs, each imparting 10 keV.

These modes were evaluated across three spatial regions:

• Region 1 (R1): A concentrated area, where the PKA energy—whether from a single atom or distributed across multiple PKAs—is localized.
• Region 2 (R2): A broader region, though some PKAs may still collide or interfere with one another.
• Region 3 (R3): A dispersed region where PKAs act almost independently, simulating isolated collisions.

This approach positions our simulations as an enhanced tool for understanding radiation-induced dynamics in graphite. Our findings indicate a particularly pronounced effect in Mode C of the dual PKAs, which consistently showed the highest number of survived FPs. The behavior seems contingent on the layered structure of graphite. When PKAs traverse different layers versus moving within a single layer, discrepancies in defect formation arise. This also can be weighed by the fact that the average TDE of graphite is around 35 eV [28]. Thus, each of the collided PKA already carries more energy than the TDE of the graphite. Thus, ultimately, more FPs can survive. Research and simulations alike have demonstrated that the bombardment by energetic clusters on various substrates can result in the creation of either craters or hillocks [51,52]. In the context of our study, exposing pristine surfaces to 40 keV cluster bombardments led to alterations in surface topography, as observed in MD snapshots. The peak intensity of the cluster beam, involving clusters of about 200–300 atoms, corresponds to average kinetic energies per atom in the range of 130–200 eV. These energy levels exceed the threshold displacement energy (TDE) for graphite, which is 35 eV.

However, the density of the energy deposited at the impact spot is relatively higher than that when it was a single PKA; thus, there is a higher chance for the local excavation of material [53]. Comparably, simulations have shown that cluster–surface collisions lead to the compression of material at the initial stage of impact and the pressure can locally rise to a GPa level [54]. Thus, more PKA would relatively lead to a stronger compression rate and, thus, a larger TDE due to compression strain and, finally, lower FPs, as shown in modes A and B as compared to mode C (Figure 3a).

Interestingly, the conventionally anticipated cascade from a single 80 keV PKA (Mode D) yielded the least defects. This can be attributed to the cascade mechanism—initial collisions might not fully expend the PKA’s energy, causing subsequent collisions to dissipate more energy as heat, reducing defect formation as the collision process itself is an instant event and energy will be dissipated at each collision until the PKA atoms stop. In Figure 4a, Modes A, B, and C exhibit a comparable peak phase diagram, characterized by an “oblong diamond-like” shape for A and B and a “diamond-like” shape for C. This configuration indicates that multiple graphite layers are affected: the broadest sections of the shape correspond to layers subjected to the most PKA energy. For instance, in Mode A, the widest layer endures the impact of six PKAs, while the topmost and bottommost layers each experience two PKAs, mirroring the distribution in Mode B. Conversely, Mode C presents a more condensed diamond-like shape configuration, with the widths of impacted layers being more uniform. This likely reduces the melted region, leading to a higher survival rate of defects when compared to Modes A and B, where a higher likelihood of reaching melting temperatures promotes greater defect recombination. Mode D, as depicted in Figure 4a, adheres to the traditional ballistic collision cascade shape, which has been reported in several works [55,56].

Regarding the regions, R1 consistently resulted in more survived defects than R2 or R3 (Figure 3). This phenomenon is attributed to thermal spikes—localized temperature surges due to energetic PKA collisions. Particularly in concentrated regions like R1, multiple collisions can elevate local temperature significantly, offering two potential outcomes. Firstly, this heightened temperature can induce atomic mobility, facilitating defect annihilation or recombination. Secondly, it may prompt increased defect generation due to enhanced atomic vibrations, leading to a net increase in defects. Figure 4b highlights the isolated nature of R3, where each PKA event operates independently without overlapping with other cascade regions. This ensures that there is no compounded effect, in contrast to the pronounced interactions observed in R1 and, to a lesser extent, in R2. Our research emphasizes that in graphite, the configuration of collision cascades is influenced not only by the energy of the PKAs, the directions in which they are displaced, and the PKA atomic type of the temperature, but also by the distribution of energy and the spatial positioning of the PKAs. Assessing multiple concurrent collision scenarios provides a more precise measure of radiation damage phases compared to examining individual PKA collisions. This insight lays the foundation for enhancing graphite designs specifically tailored for radiation environments in real conditions.

3.3. The Local-Environment Classification Scheme for Carbon Atoms Rings in Irradiated Graphite

In graphite systems, it is common to characterize atoms based on their coordination, which is generally believed to correspond to sp, sp², and sp³ hybridization, indicating two, three, and four neighboring atoms, respectively [57]. Nevertheless, observation of electron-energy-loss-spectra data for irradiated graphite has underscored the significance of categorizing threefold-coordinated atoms into four subtypes [58]. The classification of these subtypes depends on the immediate surroundings such as the coordination of second neighbors, involvement in aromatic rings, and so forth. In this work, we adopt the main four categorizations with slightly different nomenclature than described originally in this work [59]. Figure 5a depicts the atom-type composition of the four structural models, showcasing the impact of increased PKA energy and disorder. The proportion of C3-Alpha atoms (indicated by Violet color) exhibits a steady decline within the targeted irradiated damaged region. In irradiated-graphite models, the C3-Beta atoms (indicated by blue symbols) are increasing, which refer to sp² atoms belonging to a non-hexagonal ring, and which are more prevalent than C3-Alpha at 10 keV. The observed radiation damage-induced defects also include interlayer cross-links, which are known as fourfold-coordinated atoms (C4; green symbols) that are shown to be more common, especially for 40 keV and 80 keV. Unsaturated graphene edges (C2 atoms; red symbols) are responsible for the presence of threefold-coordinated atoms bonded to fourfold- or twofold-coordinated atoms [60]. These defects play a crucial role in the stability and properties of graphene-related systems [61]. The population statistics of the graphite can be viewed as a continuation of the trends seen in the other three models (i.e., C3-Alpha, and total threefold-coordinated-atom percentage decreasing with damage and the percentage of other types increasing) (Figure 5a). The total threefold-coordinated-atom percentage of graphite is somehow akin to the percentage from previous EDIP calculations of amorphous carbon [59]. The decrease in C3-Alpha and total threefold-coordinated atoms is matched to the increased irradiation damage. On the other hand, the formation of carbon rings of specific sizes is evident, as displayed in Figure 5b. Predominantly, the 7-carbon atom ring emerges as the most prevalent, with its frequency intensifying as PKA energy increases. This is followed by the 8-carbon atom ring. On the contrary, rings comprising five carbon atoms are comparatively infrequent. Notably, at elevated PKA energies, the incidence of these 5-carbon atom rings begins to diminish.

3.4. Cluster Analysis

We explored the cluster details within graphite systems, where in Figure 4b and Figure 6a, we examine the changes in the average number of interstitial clusters, while Figure 6c,d focus on vacancy clusters. One clear trend stands out: as PKA energy surges from 10 to 80 keV, there is a commensurate rise in the cluster formation in both materials. The formation of small clusters dominated by interstitial atoms, termed SIA clusters, is notably different between the two systems. In pristine graphite, the concentration of these SIA clusters, specifically in the size range of two to four atoms, is nearly half of what is evident in the hydride structure. The disparity becomes even more pronounced for clusters larger than four atoms, where the hydride structure hosts a concentration over six times that of its pristine graphite. A closer inspection reveals that pristine predominantly favors the formation of small clusters, typically comprising less than four atoms, be it vacancies or SIAs. In marked opposition, the hydride variation inclines toward the amalgamation of more sizable clusters, particularly when the PKA energy surpasses 40 keV. This highlights a key observation: the hydride structure exhibits a distinctive characteristic of generating larger cluster formations, whereas pure graphite tends to foster smaller, compact clusters. Such variances may be attributed to the inherent lattice dynamics and bonding variations existing between the two materials. The presence of hydrogen may promote the accumulation of defects due to reasons explained in Section 3.1.1, thereby facilitating the formation of larger cluster formations. Conversely, the closely knit lattice of pristine graphite might facilitate rapid defect recombination, ultimately leading to the formation of smaller clusters.

Further analysis of the data reveals a notable pattern. In pristine systems, SIA clusters are about 0.2 times more prevalent than vacancy clusters for sizes exceeding three atoms. This observation resonates with findings from other crystal studies, such as those conducted on W-Re [62] and Fe-Cr alloys [63,64]. Yet, this trend takes an unexpected turn for smaller clusters of size 2 atoms. Here, we observe that vacancy clusters outpace SIAs by a similar margin, but this reversal is specifically observed when PKA energies spike to as high as 80 keV. Such a phenomenon might be attributed to the stochastic nature of collision dynamics at elevated energies such as the focused collision sequence. The average count of vacancy clusters is relatively lower than SIA clusters, particularly for larger clusters encompassing four atoms or more. This difference can be linked to the higher formation energy of a single vacancy in graphite (7.6 eV) compared to the formation energy of an interstitial atom (6.3 eV) [65]. Specifically, the largest clusters we noticed comprised 30 SIA and 15 vacancies. Although the latter appeared more frequently, it was significantly rarer compared to other smaller-sized clusters, like groups of two or three atoms.

Upon transitioning to the hydride configuration, both the SIA and vacancy clusters exhibit comparable quantities and configurations. The inclusion of H may have an impact on lattice dynamics, thus creating an environment that facilitates the formation of both types of clusters similarly, which underscores the fact that even slight compositional alterations can have on the dynamics of defect formation in graphite.

Furthermore, we observed a distinct concave upward pattern in the formation of vacancies and SIA clusters in hydride Figure 4b and Figure 6a. When examining the rate of interstitial cluster formation, we noted that the hydrated structures showed a faster increase than the pristine graphite structures, especially at higher PKA energies.

3.5. Extra Analysis of Layered Responses to Consecutive Collision Cascades in Graphite

The same PKA was subjected to three consecutive collision events. Each of these collisions had an energy of 80 keV and spanned over 20 picoseconds for each collision (60 ps total evolution time). The aftermath of the initial collision saw a formation of 198 defects. This number increased to 364 after the second hit, but intriguingly, the second collision, in a few ps, caused an astounding rise in FPs again, tallying up to 3890 FPs, though it later stabilized to a final count of 1300 (Figure 7 and Video S6). The layered structure of graphite provides an explanation for this behavior. Each layer of graphite, while connected to the entire structure, reacts somewhat autonomously to external disruptions. Initially, the collision’s energy concentrates in the middle of the structure. But, with subsequent impacts, this energy disperses, affecting more layers but with a lessened intensity. It resembles the ripple effect in water: the initial splash is intense, but as the ripples spread out, they diminish in strength.

Post the second collision, the rapid changes observed suggest that the layers of graphite do not operate in isolation. When one-layer experiences a change, it can influence its neighboring layers, leading to a domino effect throughout the entire structure.

Interestingly, after 10 ps of the third collision, the number of FPs started to decrease. This suggests that graphite has a natural self-recovery capability. Past studies have also observed graphite’s ability to restore its structure, especially when subjected to high-temperature annealing [66,67]. However, it is crucial to note that the graphite still retains some residual stress [42], which can induce other defects if observed over an extended period. After experiencing multiple collision events, the graphite’s structure is no longer pristine. These defects can make the structure more vulnerable to future breakdowns, even though a significant portion of FPs revert to their initial configuration.

It is noteworthy that irradiated graphite shows a decrease in total threshold displacement energy. For instance, the threshold energy in pristine graphite reduces from 24 eV to approximately 5 eV in a fully disordered graphite structure [68]. This aligns with our observation of a larger volume of FPs in hydride compared to pristine graphite. The accumulation of defects alters the microstructure, diminishing its crystalline nature [69,70], and fostering cross-links between layers. Such alterations significantly influence physical attributes like thermal conductivity [71] and Young’s modulus [58]. Moreover, neutron irradiation can cause notable size changes in the graphite due to radiation damage accumulation [72].

4. Conclusions

Graphite, pivotal in many advanced technologies, exhibits vulnerability to radiation-induced modifications, especially within nuclear reactors. Employing MD simulations through the Reactive Empirical Bond Order, we delved into the atomic-level dynamics of graphite under varied irradiation scenarios, elucidating the implications of hydride inclusion. Radiation-induced defects in graphite are not merely influenced by differing PKA energy levels (10, 20, 40, and 80 keV), unique displacement trajectories, PKA species, or thermal interactions. They also hinge on the nature of the cascade patterns—whether isolated or cumulative—and collision densities. Empirically, for both graphite structures, irradiation led to an augmented defect count in the hydride compared to pristine graphite. However, both variants displayed an escalation in FPs as PKA energy surged. PKA types—either C or H—had differential impacts, with carbon inducing pronounced effects on FPs due to its greater momentum. However, hydrogen’s diffusivity varied, especially in the [0001] displacement direction, resulting in elongated diffusion times, and heightened layer penetration. Analyzing non-traditional collision modes beyond a standard single collision revealed that Mode C and the R1 region within the cascade pattern predisposed graphite to an elevated defect tally, with distinctive defect shapes such as diamond-like and oblong diamond-like morphologies in distinction to the traditional cluster-like shape as in Mode D. Our study emphasizes the importance of understanding three-coordinated graphite rings, with their nuanced atomic behaviors considering escalating PKA energies. The population statistics of graphite show a decrease in the percentage of threefold-coordinated atoms and an increase in other types of defects. The formation of carbon rings of specific sizes is also observed, with 7-carbon atom rings being the most prevalent. From a cluster perspective, SIA clusters were more prevalent in pristine graphite compared to vacancy clusters, a trend that was nearly balanced in the hydride version.