A review of displacement cascade simulations using molecular dynamics emphasizing interatomic potentials for TPBAR components

Introduction

Radiation damage in materials involves an atomic collision process initiated by high energy particles resulting in permanent defect production in the bulk or surface of the material¹. This process might include electronic excitations, electron-phonon coupling, photons, and heating of the material¹. Depending on the nature of material, crystalline materials can have bulk defects ranging from Frenkel pairs and point defect clusters² to amorphized regions³, dislocation loops^4,5 and surface defects like adatoms⁶, craters⁷, holes, bumps and ripples⁸. Materials in amorphous state on the other hand can have defects like low density regions⁹ and under coordinated regions¹⁰.

Computer simulations have been widely used to estimate and compare radiation damage in many different materials for decades¹. The production of primary damage from a single atomic collision of an energetic particle, along with the associated microstructural evolution in a material due to continuous bombardment by energetic particles, spans a wide range of length and timescales. Since material phenomena can vary across the full range of scales, which span from femtoseconds (10^-15s) for atomic collisions to years for long-term microstructural evolution, and from the atomic scale (Å) to macroscopic dimensions extending to millimeters or more, it is impossible to have one single computational method that can model material phenomena across all these scales. Accordingly, a diverse suite of computational tools is utilized to model radiation-induced phenomena based on their respective length and time scales in Fig. 1^1,11. Although mentioned in context to radiation damage, these methods are widely used for investigating other material properties and process effects apart from radiation.

**Fig. 1: Schematic representation of simulation methods by time and length scales for studying radiation effects.**

Damage can propagate in materials under three regimes based on the time scale, the primary stage, the dynamic recombination stage and the steady state^5,12. The primary damage stage (0.1–1 ps) or the ballistic stage is characterized by rapid atomic collision and temperature spike caused by non-equilibrium dynamics^13,14. The dynamic recombination stage (1–10 ps) is characterized by the formation and annihilation of point defects and their clusters. The steady state is reached by the surviving defects that are not annihilated in the recombination stage¹⁴. The surviving defects continue to evolve over a larger timescale of years resulting in degradation of material properties.

The different simulation techniques in Fig. 1 can be linked to the physics of different stages of radiation damage. Radiation damage in the bulk material is a result of energy transfer from a passing high energy particle such as an ion or a neutron. The penetrating particle interacts and transfers its high energy to an atom in its path. If the transferred energy is higher than the threshold displacement energy (TDE), then the energized atom itself becomes the primary knock-on atom (PKA) that travels ballistically and collides with other atoms in the lattice producing secondary knock-on atoms. The sequence of repeated collisions between atoms results in a collision cascade that results in primary damage in the material. The initial stage of cascade collisions occurs on a ~10 femtoseconds, which is at least an order of magnitude faster than the timescale of lattice vibrations (300 fs) which makes it a non-thermodynamic process resulting in a thermal spike¹⁴. The initial stage of the radiation damage can be simulated by molecular dynamics (MD) or binary collision approximations (BCA) as they well encompass the required timescale. For the generated defects, DFT simulations can be useful in calculating their migration properties (typically by using the nudged elastic band method (NEB))¹. Subsequently, Kinetic Monte Carlo (KMC) calculations can be used to calculate the diffusion rates. KMC and discrete dislocation dynamics (DDD) can be used further to simulate dislocation mobility at the mesoscale¹⁴.

Tritium, a critical hydrogen isotope, is produced through the neutron irradiation of Tritium Producing Burnable Absorber Rods (TPBARs), which are specifically engineered for both the production and capture of tritium^15,16. A schematic of a TPBAR is depicted in Fig. 2¹⁷. Each TPBAR contains pellets made of lithium aluminate, enriched with ⁶Li, which generates tritium when exposed to neutron radiation. The tritium produced is subsequently absorbed by a surrounding Zircaloy-4 getter tube that surrounds the LiAlO₂ pellet¹⁸.

All these components are housed within a stainless steel cladding, to inhibit the diffusion of tritium outward into the reactor coolant. The inner surface of the cladding is coated with Ni. This review aims to examine the MD methods for displacement cascades reported in the literature, pinpointing critical elements of simulation techniques, and describing the current state of the art. The focus is specifically on materials used in TPBARs and providing an in-depth analysis of cascade simulations relevant to these materials¹⁵. Table 1 presents a list of materials relevant to TPBAR components investigated in numerous studies, as compiled from the literature.

Table 1 Summary of cascade simulations using molecular dynamics across different materials relevant to the TPBAR and some other nuclear materials of interest

Full size table

Dependence of peak and surviving defects on interatomic potentials

When considering peak damage versus the surviving damage, the nature of the material and PKA energy play a prime role. Most works in the literature are based on single crystal evaluations and use Frenkel pairs to quantify damage. Frenkel pairs are a fundamental type of point defect that consists of a vacancy, where an atom is missing from its lattice site, and a corresponding interstitial, where the displaced atom occupies a non-lattice position¹. Beland et al.¹⁹ showed that in Ni the peak in the number of defects at 40 keV PKA is around 5000 while that in NiCo and NiFe is higher up to 10,000. But the number of surviving defects is lower in alloys than in pure Ni, in the order of a 300–500 after 15 ps. The lower number of stable defects in the alloys than in the pure Ni was attributed to the harder interactions between the Fe-Fe, Co-Co, Ni-Fe and Ni-Co pairs than the Ni-Ni interactions. The peak damage in an equiatomic quinary alloy FeNiCrCoCu was found to be around 6000 FPs at 40 keV PKA with fewer than 100 surviving FPs²⁰. HCP Zr²¹ showed a peak damage of 6000 FPs at 80 keV and a surviving damage of ~100 FPs after 35 ps. CoCrFeMnNi HEA showed a peak damage of ~ 25000 defects at 10 keV at 773 K as compared to 6000 for Ni at the same PKA energy and temperature²² but this work used a MEAM potential as compared to the EAM potential used by Beland et al.¹⁹. The plots of the numbers of Frenkel pairs versus time are compared in Fig. 3, from 3 different works investigating pure Ni using different potentials. Comparing the black and the green lines, representing 20 keV PKA cascades at 100 K and 300 K respectively, higher temperature should ideally cause a higher peak damage due to thermal activation. However the black line²³ at 100 K shows a higher peak damage than the green line¹⁹ at 300 K. The number of surviving defects on the other hand should be lower at higher temperatures^24,25 but the investigation in ref. ²² (blue line) representing 10 keV at 773 K contradicts that in ref. ¹⁹ (red line) representing 10 keV at 300 K with a slightly higher number of surviving defects at higher temperature.

**Fig. 3: The variation in profiles of the Frenkel pairs plotted as a function of time for pure Ni as studied by Voskoboinikov et al.**

This emphasizes the criticality of the potentials chosen for cascade simulations, where there can be an order of magnitude mismatch in the defects produced in the same material if the potential is switched from one to another. A NiCoCrFe alloy shows a peak damage of ~5000 FPs and a surviving damage of <100 FPs which closely matches that in Ni at 20 keV²⁶. However, there is a 1% higher recombination in HEA than in pure metals²⁶. The effect of temperature on cascade damage in Ni is elaborated in ref. ²³ where high temperatures of up to 1200 K in Ni produced lower surviving damage of around 30 FPs while at lower temperatures the FPs were as high as 70.

The effect of using electronic stopping or the two-temperature model in MD can swing results by a factor of 5–6 in some metals and alloys²⁷. For example in NiFe alloy at 50 keV PKA, shows ~100,000 peak FPs with an average of 20,000 surviving defects in classical MD, but the number peak FPs falls by almost a factor of 6 if electronic stopping is employed in the simulations²⁷. On the other hand, in NiFeCr there is no significant effect of electronic stopping²⁷. Electronic effects in damage cascade simulations should be included when studying materials with a high atomic number (Z) or when simulating high-energy collision cascades²⁸. In such cases, a significant portion of the kinetic energy of the primary knock-on atoms (PKAs) is transferred to the electronic subsystem, which affects the dynamics of the cascade and the resulting defect structures. The inclusion of electronic stopping, typically modeled as a friction term in the equations of motion, modifies the atomic trajectories by dissipating energy into the electronic subsystem. This dissipation depends on both the atomic mass and Z-value of the material²⁸. For example in a SiC system, C atoms with a lower mass than Si atoms, will have a higher velocity for a given kinetic energy and will thus have more electronic stopping that Si atoms with a higher atomic mass²⁹.

While considering materials like Si, the peak defects produced by a 5 keV PKA is 600 Frenkel pairs, out of which fewer than 100 survive eventually³⁰. However, for a 10 keV PKA the peak was higher at 6000 defects with about 1050 surviving defects as shown in ref. ³¹. This discrepancy of an order of magnitude is again likely to originate from the different potentials employed in the two works, where the former uses a Stillinger-Weber potential while the latter uses a Tersoff potential. Samolyuk et al.³² drew an excellent comparison of Tersoff and Gao-Weber (GW) potentials in cascade simulations of SiC. Specifically, they found that the ratio of peak number of defects to the surving defects is considerably higher for GW potential as compared to Tersoff potential. This was because the barrier for carbon interstitial-vacancy recombination is predicted to be higher by the Tersoff potential than by the GW potential³².

In ceramics like LiAlO₂, the numbers of peak and surviving defects at 15 keV were found to be in the range of 200 and 120 respectively¹³. While in the secondary phase of LiAlO₂, formulated as LiAl₅O₈ the numbers of peak and surviving defects at 15 keV were found to be ~60 and ~45 respectively. These ceramics are a critical component in tritium-producing burnable absorber rods (TPBARs) where they are exposed to neutron irradiation to produce tritium³³. A challenge faced in simulating such a system are the core-shell model used for oxygen which poses the risk of core-shell separation at high PKA energies. Another work³⁴ used a core-shell model for oxygen in Li₂TiO₃ and obtained a peak defect count of upto 1500 and surviving defects up to 160 for a 2 keV PKA.

Diamond is an important component in radiation detectors for particle accelerators^35,36 and has a high tolerance to radiation damage as shown in ref. ³⁷. Only 22 FPs are produced at the peak damage at 2 keV and around 10 surviving FPs³⁷. Diamond showed a smaller PKA range and cascade length compared to graphite³⁸ indicating a higher tolerance to radiation damage than graphite.

Criteria for selecting the optimal interatomic potential for cascade simulations

Interatomic potentials determine the energy change during atom collisions and subsequent atomic displacements. They affect the diffusion and migration of defects created during the cascade. These behaviors are essential for comprehending the long-term stability and evolution of materials under radiation. Different potentials can result in varying spatial distributions of damage. A notable work by Terentyev et al.³⁹ compares four different interatomic potentials for Fe— the long-range EAM potential denoted by AMS⁴⁰, a different version of the EAM denoted by WOL⁴¹, another long-range EAM denoted by CWP⁴², and the short-range Finnis-Sinclair type potential denoted by ABC⁴³—focusing on the number of defects produced, cascade density, and the behavior during and after the cascade. The main observations can be summarized as follows: in terms of peak defects and cascade density, WOL produces significantly fewer defects at peak time with a lower cascade density due to a larger average cascade volume. AMS generates the highest number of defects at peak time (2–5 times more than WOL) and has the highest cascade density. ABC and CWP fall in between, with ABC closer to AMS and CWP closer to WOL. Subcascades are observed mainly above 20 keV cascade energy, aligning with Stoller’s findings⁴⁴. The CWP potential shows a stronger tendency toward subcascade formation compared to the other potentials. However, WOL produces cascades that are too diffuse to form subcascades. Despite differences in peak defect numbers, all potentials exhibit a similar efficiency (around 0.3 ± 0.1) in terms of surviving Frenkel pairs (FPs), which is consistent with previous studies and experimental data. An exception is WOL, which has slightly higher efficiency despite having the fewest defects at peak. The stiffness-to-range (S/R) ratio at ~ 30 eV correlates with defect production. Higher S/R ratios, as seen with WOL, result in more dilute cascades and potentially higher numbers of surviving defects due to reduced recombination. In contrast, softer and longer-ranged potentials, like AMS, favor dense cascades and more efficient defect recombination. A striking conclusion in this work was that there was no correlation between the frenkel pairs produced and the TDEs.

Another work by Gao et al.⁴⁵ investigated the influence of different interatomic potentials on displacement cascade simulations in α-Fe. One of the key findings is that the TDEs did not correlate with the number of Frenkel pairs produced in the cascades. The work compared the performance of three magnetic potentials (MP) from refs. ^46,47 and a Mendelev-type potential⁴⁰, with the magnetic potentials being hardened to make them suitable for cascade simulations. The results showed that the peak time, maximum number of defects, cascade volume, and cascade density with the magnetic potentials were smaller than those observed with the Mendelev-type potential. Interestingly, there was no significant difference, within statistical uncertainty, in the defect production efficiency between the Mendelev-type potential and the second magnetic potential at the same cascade energy. However, the defect production efficiency was remarkably smaller with the first and third magnetic potentials. The study also revealed that the self-interstitial atom (SIA) clustered fractions were smaller when using the magnetic potentials compared to the Mendelev-type potential, particularly at higher energies. This behavior was attributed to the larger interstitial formation energies produced by the magnetic potentials, leading to fewer defect clusters at higher cascade energies. The defect clustered fractions, which are critical input data for radiation damage accumulation models, could significantly impact predictions of microstructural evolution under irradiation. As the cascade energy increased, so did the peak time, maximum number of defects, cascade volume, and cascade density, but these quantities were generally smaller with magnetic potentials compared to the Mendelev-type potential. The study also observed that larger interstitial formation energies from the magnetic potentials resulted in fewer defects at the peak time and smaller cascade volumes and peak times. While cascade density increased with energy, it appeared to be independent of the stiffness-to-range (S/R) ratio. Interestingly, there was no direct correlation between the TDEs and the number of defects surviving the cascades. The hardening procedure of the magnetic potentials and the Mendelev-type potential may have contributed to differences in the number of Frenkel pairs produced during displacement cascades. The study concluded that the SIA clustered fractions were smaller for the magnetic potentials due to the higher interstitial formation energies, while the vacancy clustered fractions were largely similar across potentials, except at higher energies where the Mendelev-type potential produced larger fractions. Overall, these differences in defect clustered fractions may influence microstructural predictions in radiation damage accumulation models.

Sand et al.⁴⁸ conducted an evaluation of five interatomic potentials for tungsten, including the Finnis-Sinclair potential⁴⁹ modified by Ackland and Thetford⁵⁰, and by Zhong⁵¹ (ATZN); the potential by Derlet et al.⁵², modified by Björkas et al.⁵³ (DNDBN); the potential by Justin et al. (JW)⁵⁴; and a new potential developed by Marincia et al.⁵⁵, stiffened by Sand et al.⁴⁸ through soft (MS_soft) and hard (MS_hard) interpolation. The study found that all five potentials produced similar predictions for vacancy formation energy, with values close to the experimental result of 3.7 eV. However, significant variation was observed in the calculated TDEs across different crystallographic directions. For the <100> direction, the TDEs ranged from 31 to 63 eV, with the experimental value being 40 eV. In the <110> direction, TDE predictions ranged from 51 to 103 eV, while the experimental value was 70 eV. Similarly, for the <111> direction, the predicted TDEs spanned from 41 to 89 eV, compared to the experimental value of 44 eV. Despite these variations, the potentials that were otherwise similar in their characteristics exhibited good correlation in their predictions of defect numbers for cascades initiated by primary knock-on atoms (PKAs) with energies below 50 keV, based on their TDEs. Interestingly, the study also found that there was no correlation between the final number of defects produced in the cascades and the defect formation energies predicted by each potential. This suggests that, while TDEs play an important role in predicting defect numbers at intermediate energies, the relationship between defect formation energies and the final damage state in higher-energy cascades may be more complex.

Therefore, when selecting an interatomic potential for cascade simulations, the TDE is a critical parameter to consider. TDE defines the minimum energy required to displace an atom from its lattice site and initiate a defect. A potential with accurate TDEs across different crystallographic directions is crucial because it directly influences the initiation of displacement cascades. If the potential underestimates TDEs, it may predict an unrealistic number of defects at low energies, while overestimating TDEs could lead to underprediction of damage.

When multiple machine-learning potentials are available for a system, such as GAP and tabGAP for tungsten⁵⁶, it’s essential to balance computational efficiency and accuracy. One effective metric for gauging accuracy, especially in the context of radiation damage simulations, is the TDE. TDE directly relates to the primary damage stage and offers a reliable criterion to assess the potential’s suitability for cascade simulations. Therefore, choosing a potential should involve evaluating both its computational cost and performance in critical quantities like TDE, ensuring an optimal balance between speed and accuracy.

Impact of short-range repulsive potential splining

In radiation damage simulations, modifying interatomic potentials to accurately describe repulsive interactions at short distances is essential for capturing both the ballistic and thermal phases of a displacement cascade. A common approach involves using the universal ZBL potential to represent repulsive interactions at interatomic distances below 1 Å. The transition between the ZBL potential and the equilibrium potential is typically achieved via a smooth interpolation function in the intermediate repulsive range, which can be tuned to reproduce accurate TDE. This combination allows for modeling highly repulsive interactions during collision events, while also maintaining accurate near-equilibrium behavior. Sand et al.⁴⁸ employed a fifth-degree polynomial to spline the ZBL potential and created two interpolations for a tungsten potential by Marinica⁵⁵ as shown in Fig. 4 (a). The M-S_h (stiffness S = -140.5 eV/Å, R = 1.5 Å, S/R = 93.7 eV/Å²) produced relatively lower atomic displacements by the end of the ballistic phase (~300 fs), with additional atomic mixing continuing during the thermal phase, driven by the heat spike. Conversely, the M-S_s potential (S = -552.0 eV/Å, R = 1.1 Å, S/R = 501.8 eV/Å²) resulted in longer ballistic recoils and less energy transfer to surrounding atoms, leading to a more diffuse cascade without significant heat spike diffusion as shown in Fig. 4 (b). The M-S_s potential caused more surviving defects due to reduced in-cascade recombination, demonstrating how potential stiffness directly influences cascade evolution and defect survival.

**Fig. 4: Comparison of potential energy and displacement dynamics using various interatomic potentials.**

Byggmästar et al.⁵⁷ conducted a comparative analysis of two interatomic potentials for iron, focusing on the impact of the ZBL spline modifications on cascade simulations. They evaluated the Marinica^58,59 potential (M07) and its modified version (M07-B) using quasi-static drag (QSD) simulations, where an atom is stepwise moved through the crystal lattice, testing many-body interactions and comparing the results to ab initio data⁶⁰. Through this process, they fine-tuned the potential to achieve both satisfactory TDEs and a good agreement with DFT energy curves. The M07 potential (S = 101 eV/Å, R = 1.40 Å, S/R = 72.3 eV/Å²) was adjusted to produce the stiffer M07-B potential (S = 143.6 eV/Å, R = 1.30 Å, S/R = 110.5 eV/Å²). This modification in the stiffness of the repulsive part of the potential significantly affected the cascade characteristics. The M07-B potential generated a higher number of Frenkel pairs at the same PKA energy compared to M07, with peak damage occurring earlier and the relaxation time extending longer. M07-B also produced larger vacancy clusters than M07, which was attributed to differences in the dynamics of the early cascade stages, controlled by the repulsive part of the potential. The stiffer spline in M07-B led to more energetic recoils and less recombination during the ballistic phase, resulting in higher defect numbers and more extensive vacancy clustering. This study highlighted how tuning the ZBL spline can influence defect production and cascade evolution, directly linking the number of Frenkel pairs to the average TDE of the system.

The works of Sand et al.⁴⁸ and Byggmästar et al.⁵⁷ highlight the critical role of splining in modifying interatomic potentials that directly influences cascade dynamics. In both tungsten and iron, adjustments to the stiffness and range of the repulsive spline impact the ballistic phase of displacement cascades. For stiffer potentials, such as the M07-B in iron and the hard-splined potential in tungsten, cascades result in more pronounced vacancy clustering and higher Frenkel pair production due to stronger recoils and less energy dissipation. Softer splines, like M07 and the softer tungsten potential, produce more diffuse cascades, leading to lower defect clustering and greater recombination of defects during the thermal spike. These studies demonstrate that the spline’s stiffness controls the recoil dynamics and the balance between energy dissipation and defect production.

Advent of machine learning potentials and differences from classical potentials

In recent years, there has been a surge in the application of machine-learning (ML) techniques to develop interatomic potentials. These methods encompass various architectures and descriptors designed to effectively capture the atomic environments effectively^61,62. Established frameworks cover neural networks^63,64, Gaussian process regression and kernel methods^65,66,67 and linear regression^68,69, each offering unique advantages. Despite being relatively new, the field has rapidly matured, enabling the routine training of high-quality machine-learning potentials for various materials and molecules^{70,71,72,73,74}. This evolution underscores the growing reliability and adaptability of machine-learning methodologies in simulating atomic interactions, heralding significant strides in materials science and computational chemistry.

The most common interatomic potential used for simulating displacement cascade simulations in metals and alloys is the EAM potential⁷⁵ which is formulated as

$${E}={{F}_{!a}}mathop{sum}limits_{{j}ne {i}}{rho}_{i}({R}_{i}.{j})+frac{1}{2}mathop{sum}limits_{{j}ne{i}}{{{varnothing}}}_{alpha ,beta }({R}_{i}.{j})$$

(1)

where E is the total energy of the atomistic system, R_i,j represents the distance between jth and ith atoms, ρ is the electron density, F represents the embedding energy of the ith atom, ∅ is the short range pairwise potential energy, and α and β being the elements of atoms i and j. The core concept of the EAM potential lies in the notion of embedding energy. Each atom in the metal lattice experiences an effective potential, or embedding energy, due to its interaction with the surrounding electron density. This embedding energy is a function of the local electron density at the atom’s position. In addition to the embedding energy, there is also a pairwise potential that accounts for the interaction between pairs of atoms. This potential captures the repulsive and attractive forces between neighboring atoms and typically follows a functional form determined empirically or theoretically.

Machine learning potentials (MLPs) can be broadly categorized into six distinct types based on the underlying machine learning model employed in their construction⁷⁶. The first category includes artificial neural networks (ANN) and deep neural networks (DNN), as originally proposed by Behler⁶³. These networks have led to the development of advanced potentials like the deep tensor neural network (DTNN)⁷⁷, hierarchical interacting particle neural network (HIP-NN)⁷⁸, continuous filter convolutional layers (SchNet)⁷⁹, and DeepPotential⁸⁰, all of which use neural networks to model complex atomic interactions with high accuracy. The second category is graph neural networks (GNN)^81,82, which utilize graph structures where atoms are represented as nodes and bonds as edges, capturing atomic interactions in a graph convolutional framework. Notable GNN-based potentials include GNNFF⁸³, MDGNN⁸⁴, TeaNet⁸⁵, GEMNet⁸⁶, HermNet⁸⁷, MACE⁸⁸, Allegro⁸⁹, and CHGNet⁹⁰. A third category, equivalent graph neural networks (EGNN)⁹¹, extends the GNN framework by incorporating symmetries such as rotation, translation, and reflection invariance. This approach allows EGNN models, such as NequIP⁹² and UNET⁹³ to effectively scale to higher dimensions while preserving the physical symmetries of the system. The fourth category includes polynomial potentials like the Moment Tensor Potential (MTP)⁶⁹ and atomic cluster expansion (aPIP)⁹⁴, which rely on polynomial fitting to represent interatomic forces. These models offer a balance between accuracy and simplicity, making them efficient for large-scale simulations. The Gaussian Approximation Potential (GAP)⁶⁵ represents the fifth category, which bridges the gap between traditional empirical potentials and first-principles methods by incorporating Gaussian process regression. Examples like FCHL19⁹⁵ and GAP are highly flexible and offer a detailed representation of atomic interactions, often incorporating quantum mechanical effects. Finally, the sixth category comprises models based on equivariant transformer architectures, such as TorchMD-NET⁹⁶. These models balance accuracy and computational efficiency, making them particularly useful for molecular dynamics simulations where fast yet reliable predictions are essential.

Datasets form the cornerstone of machine learning, offering the training information that algorithms use to learn and generate outputs. Large datasets are crucial for developing robust models capable of generalizing well across various scenarios, minimizing biases, and enhancing prediction accuracy. In materials science, gathering a wide range of datasets is essential to understanding material properties, behaviors, and interactions, which helps predict performance and discover novel materials with specific characteristics. To support MLPs, a significant method involves leveraging open-source materials databases. Initiatives like materials genomics have led to the establishment of large-scale databases worldwide, including the ICSD⁹⁷, Materials Project⁹⁸, Aflow⁹⁹, Materials Cloud¹⁰⁰, NOMAD¹⁰¹, ALKEMIE MatterDB¹⁰² and MatNavi¹⁰³. These databases provide extensive DFT-calculated data, aiding in the training of ML models without the need for time-consuming calculations.

A typical MLP however is quite different from traditional pairwise potentials. MLPs employ a regression model to predict the potential energy of a system of atoms based on their atomic configurations. The specific equation for an MLP can vary depending on the chosen machine learning algorithm, the representation of atomic configurations, and the complexity of the model. However, a general form of an MLP equation can be expressed as follows:

$$Eleft(Rright)=f(R,theta )$$

(2)

Here, (Eleft(Rright)), represents the potential energy of the system, which is a function of the atomic positions represented by the vector R. The function f is the regression model, which takes the atomic configurations R as input and outputs the corresponding potential energy. The parameter(s) (theta) represents the weights and biases of the regression model, which are learned from training data. The choice of f and the specific form of the regression model can vary depending on the nature of the material and descriptors. The atomic configurations R may be represented using various descriptors, such as atomic positions, atomic distances, radial distribution functions, or other structural features. In summary, a typical MLP equation involves a regression model f that maps atomic configurations to potential energy, with parameters (theta) learned from training data. The specific form of and the representation of atomic configurations can vary depending on the chosen machine learning algorithm and the characteristics of the system being studied.

Machine learning potentials for displacement cascade simulations

Machine learning potentials are intricate mathematical models designed to capture potential energy surfaces (PES) by using quantum mechanical data¹⁰⁴. The development of such a potential involves three crucial steps: The first is curating a database of representative structures and associating them with reference values such as energies, forces, and stresses obtained through accurate methods such as DFT and beyond¹⁰⁵. The second step involves building descriptors analogous to feature vectors in traditional machine learning. These descriptors need to exhibit invariance to atomic transformations (translations, rotations and permutations)¹⁰⁶. The final step involves training a flexible, highly parameterized function to map descriptors to reference data. Various fitting approaches are employed, including artificial neural networks (such as Behler–Parrinello-type networks^63,64,107, DeepMD¹⁰⁸, SchNet¹⁰⁹), kernel-based methods (e.g., Gaussian Approximation Potentials⁶⁵), and linear models (e.g., Spectral Neighbor Analysis Potential (SNAP)⁶⁸, Atomic Cluster Expansion (ACE)¹¹⁰). The selection of the fitting methodology, along with associated hyperparameters, significantly influences the potential’s accuracy, generalizability, and computational efficiency. Machine learning interatomic potentials provide a significant advantage over conventional empirical potentials in MD simulations across various material systems. Traditionally, interatomic potentials have been derived from empirical methods, such as the EAM, empirical N-body potentials, and Tersoff potentials, which offer computational efficiency but often lack the accuracy and transferability needed for complex systems. These methods rely on predefined functional forms and parameters fitted to limited datasets, which can restrict their predictive capabilities, especially in diverse environments. ML methods, which can learn directly from large datasets and capture intricate relationships in atomic interactions, offer a compelling alternative to both empirical and first-principles methods. ML-based potentials improve accuracy by leveraging extensive training data and offer greater transferability across a wide range of materials and conditions. This paradigm shift holds promise for advancing the precision and scalability of materials simulations^76,111.

With regards to developing ML potentials for cascade simulations, major contributions to the community have been made by the University of Helsinki, Finland^112,113,114. Numerous potentials have been created for metals such as Mo, Ta, V, Nb, and W,^{50,52,115,116} effectively replicating a range of properties with satisfactory precision while remaining computationally efficient, they have exhibited limitations in accurately reproducing certain quantities. These limitations include surface energies, dislocation energies, as well as vacancy and self-interstitial cluster energies¹¹⁷. As a response to these challenges, Gaussian approximation potentials (GAP) were specifically crafted for metals such as V, Nb, Mo, Ta, and W¹¹⁷. These potentials have demonstrated notable accuracy for targeted properties, encompassing fundamental elastic and thermal characteristics, defect and surface energetics, and liquid-phase properties. They also exhibit robust transferability to properties beyond those directly addressed by the structures in the training database. For instance, they accurately predict the formation energies of self-interstitial clusters and the melting behavior under extreme pressures, both of which are critical in simulations pertaining to radiation damage¹¹⁷. Recent reports have highlighted refractory high entropy alloys to be good candidates for nuclear applications. However the challenges arise due to random distribution of large number of constituent elements (>5) that result in corrugated distribution of energy potentials leading to a complex energy landscape. This distorted lattice structure¹¹⁸ may help in raising the defect formation and migration energy^119,120 or TDE and improve the radiation damage tolerance of HEAs^121,122. Fan et al.¹²³ developed unified general MLP for 16 elemental metals and their alloys, demonstrating an efficient method to represent the vast chemical space without generating training data for every possible combination. This work highlights the application of MLPs in exploring the mechanical behavior of complex materials, such as plasticity and primary radiation damage in MoTaVW refractory high-entropy alloys. Their approach, utilizing distinct neural networks for each atomic species and multiple loss functions for parameter optimization, is both scalable and adaptable, suggesting a pathway for future models capable of covering the entire periodic table.

Recently, the suitability of the GAP potential for tungsten (W) was specifically investigated in the context of cascade simulations¹¹³. In particular, the authors showed that the machine-learning potential (GAP) trained on a moderately sized dataset demonstrates the capability to accurately represent various properties of tungsten akin to the accuracy of DFT precision. Although the potential has a general nature, the researchers particularly focus on replicating properties pertinent to radiation damage. The adaptable nature of the machine-learning framework enables the potential to effectively describe properties that have historically posed challenges for analytical potentials, including the stability of defect clusters and diverse surface characteristics. Consequently, this ML potential proves valuable for extracting more precise data from classical molecular dynamics simulations, particularly in the context of radiation damage in tungsten relevant to fusion applications. By doing so, it facilitates the resolution of previously ambiguous disparities observed in the results obtained using different existing potentials¹²⁴. The authors also tested the potential on a range of tungsten properties such as surface energies, bulk properties like twinning energies and energy-volume curves, as well as phonon dispersion and self-interstitial atoms and cluster energies as a function of size as shown in Fig. 5¹¹³. Encouragingly, they found good agreement with values obtained from DFT simulations for these properties. However, it’s essential to note that the computational cost associated with this GAP implementation is considerably higher compared to traditional analytical potentials, by ~2–3 orders of magnitude. This increased computational demand poses a challenge in acquiring extensive statistics on the primary damage. Nonetheless, recent advancements in optimizing the smooth overlap of atomic positions (SOAP) kernel¹²⁵ have shown promising enhancements in computational efficiency without compromising accuracy.

**Fig. 5: Formation energies of self-interstitial clusters obtained from GAP and compared with DFT.**

This GAP potential was further advanced in a follow-up work to produce a faster performance version called the tabulated GAP (tabGAP)¹²⁶. The tabGAP achieves a computational speed that is two orders of magnitude faster than the original GAP, while maintaining a similar performance in terms of the number of surviving defects and sizes of defect clusters. TabGAP operates by pre-computing GAP energy predictions and mapping them onto low-dimensional grids for rapid access during simulations. This tabulation allows for the use of conventional spline interpolation methods, significantly accelerating computation.

Koskenniemi et al.⁵⁶ used tabGAP to simulate primary radiation damage in 50-50 W-Mo alloys and pure W using MD. The results indicate that W-Mo alloys exhibit a defect survival rate comparable to that of pure W and show more effective defect recombination during the initial cascade phase. In certain scenarios, W-Mo alloys also show complete defect recombination post-cascade cooling, which is not observed in pure W.

Despite the increased speed, tabGAP is still two orders of magnitude slower than traditional EAM potentials. An optimized version of tabGAP, with improved code and adjusted cut-off radii, achieves speeds three orders of magnitude faster than GAP and only one order slower than EAM potentials. The defect clustering performance of tabGAP closely mirrors that of GAP, within standard error margins, although minor differences in specific cluster sizes were noted⁵⁶. This balance of computational efficiency and accuracy in simulating radiation effects makes tabGAP a valuable tool in the field of MD simulations.

Fission products inclusion in cascade simulations

It’s important to emphasize that irradiation processes in TPBARs generate tritium and helium as byproducts, both of which contribute significantly to radiation damage in the pellets and surrounding materials. Similar to fission gas behavior in traditional fuels like UO₂, these gas byproducts in TPBARs can result in bubble formation, pressure build-up, and subsequent material degradation due to internal stresses. Fission gas evolution simulations using molecular dynamics are computationally challenging in the sense that the interatomic potentials defining the fission gas and fuel/component interactions are very limited in the literature. While gas interaction in lithium aluminate has not been extensively studied, prior work has examined gas bubbles in fuels and materials. One such combination is the He and Xe in UO₂ that has been investigated in the several papers^{127,128,129,130}. MD simulations were used in a¹²⁷ study on helium and xenon gas bubble behavior in UO₂. The work utilized the Basak potential¹³¹ for U-U and U-O interactions, the Geng potential¹³² for Xe-UO₂ interactions and the Grimes potential¹³³ for the He-UO₂ interactions. Although this work¹²⁷ did not simulate displacement cascades in UO₂ fuel, it highlights the pressure release mechanisms of He and Xe bubbles. In the case of helium bubbles, pressure relief occurs through the creation of defects on the oxygen sublattice. This facilitates helium atoms to diffuse into the adjacent bulk UO₂ from the bubbles, forming a diffuse interface. Conversely, xenon bubbles release pressure by displacing uranium atoms, leading to U aggregation around the xenon bubble as a pressure release mechanism. An earlier work by the same group¹²⁸ used the same set of interatomic potentials to model Xe diffusion in UO₂. Under nonequilibrium conditions, xenon atoms occupy octahedral interstitial sites instead of the stable uranium substitutional site. Utilizing basin-constrained molecular dynamics the work identified a transient migration mechanism involving xenon and two oxygen atoms¹²⁸.

Cascade simulations on UO₂ in the presence of the fission gas Xe were a one kind of study performed by ref. ¹³⁰. However, this work used a different potential than the work described in the previous paragraph. Here, U-O interactions were described by the Morelon potential¹³⁴ similar in format to the Buckingham potentials, and the Xe-Xe interactions were defined by the resolution of the Geng potential¹³². The work investigated the Xe bubble resolution into the UO₂ matrix. The resolution process involved temporary trapping of fission gas by intragranular bubbles in nuclear fuel, which was subsequently brought back into atomic solution through interactions with fission fragments or fast neutrons. The simulations revealed that low-energy PKAs were ineffective in the resolution process. High-energy interactions destroyed smaller bubbles entirely and maintained a quasi-constant number of gas atoms in resolution when interacting with larger bubbles.

Another work investigating He resolution from bubbles in iron was performed using MD by ref. ¹³⁵. A mix of interatomic potentials was employed, Fe-He was defined by the Oak Ridge potential^136,137, the iron matrix was described by the Ackland potential⁴³ and He-He interactions were defined by the Aziz potential¹³⁸. The work used MD simulations to explore variables such as the irradiation temperature (100 and 600 K), cascade energy (5 and 20 keV), bubble radius (0.5 and 1.0 nm), and the ratio of helium to vacancies within the bubble (0.25, 0.5, and 1.0). Consistent patterns emerged, indicating that larger bubbles and those with higher helium-to-vacancy ratios experienced more significant helium displacement through ballistic resolutioning. The efficacy of resolutioning decreased at higher temperatures (600 K) compared to lower temperatures (100 K) and with higher energy cascades (20 keV) compared to lower energy ones (5 keV). Overall, the findings suggest a fair extent of helium removal via ballistic resolutioning, with potential implications for fusion irradiation scenarios characterized by elevated helium levels resulting from transmutation processes.

Recommendations for new interatomic potentials for TPBAR components

An important intent of this review paper is to provide recommendations for the appropriate potentials available to model each material of the TPBAR components using displacement cascades in MD. Each material will be discussed further, and the past relevant works will be highlighted to show how they inform current best practices in interatomic potential selection. This analysis is critical for enhancing the accuracy and reliability of simulations intended to optimize TPBAR performance and safety.

SS 316 for cladding

In the early 2000s, a few works attempted to study steels by using just Fe-Cr based compositions due to lack of interatomic potentials that could define interactions between all elements of steel, namely Fe, Cr, Ni, Mo, C, Mn, P, S, Si, N. Although Fe and Cr have the highest wt. % in the composition, the presence of the other elements in steel significantly tunes its properties. Displacement cascades were simulated in Fe-10%Cr alloys up to 15 keV¹³⁹ and then later up to 50 keV¹⁴⁰. In the later study¹³⁹ it was found that the addition of 10% chromium does not significantly alter the initial collisional phase of the cascade, compared to pure Fe. However, it does influence the post-collisional stage, specifically affecting the redistribution of dumbbell species. A notable observation was that most interstitial atoms are chromium, present in a much higher proportion than its overall concentration in the alloy. This disproportionate presence of chromium appears to inhibit the recombination process, consequently leading to a slight increase in defect production in the alloy when compared to pure Fe as shown in Fig. 6. No marked differences were observed in the proportion of defects that cluster between the two materials. Nevertheless, clusters in the Fe–10%Cr alloy predominantly contain chromium atoms, which seem to stabilize these clusters. This stabilization significantly reduces the mobility of interstitial loops in the alloy compared to those in pure Fe.

**Fig. 6: Distribution of Fe-Fe, Fe-Cr, and Cr-Cr dumbbells in Fe-10%Cr during cascade simulations.**

The latter study¹⁴⁰ investigated 50 keV PKA cascades in Fe-10%Cr and found that the concentration of chromium (Cr) in interstitial defects increases predominantly in the post-collisional and post-relaxation stages. This buildup occurs as single interstitial atoms move until they encounter and become immobilized by Cr atoms. Simultaneously, clusters are likely to nucleate and grow, and coalesce in areas with elevated local concentrations of Cr, where their movement is significantly hindered. This sequence of events highlights the effect of Cr on defect dynamics and defect cluster stabilization within the alloy under irradiation conditions.

In 2013, Bonny, Castin, and Terentyev¹⁴¹ developed a potential Fe-Ni-Cr that allowed the inclusion of Ni in the existing Fe-Cr composition to represent the SS-316L composition. Using this potential Collette et al.¹⁴² performed MD simulations to examine the influence of pre-existing dislocations on radiation-induced defect formation. In their simulations, edge dislocations were introduced into the crystal structure by forming planar defects along typical FCC slip systems, replicating the dislocation tangles seen in experimental TEM studies of AM SS-316 L samples. PKAs with energies of 5, 10, and 15 keV were introduced at the center of these tangles, resulting in cascade volumes that consistently intersected with the dislocation cores. Findings from the simulations suggested that areas dense with dislocations tend to lower the likelihood of surviving point defects clustering together. This reduction is attributed to the defects being absorbed by the dislocation cores, which then break down into partial dislocations.

In the context of TPBARs, the primary concern in simulations typically revolves around radiation damage. However, the permeation of tritium through the cladding is of significant importance. Consequently, understanding tritium transport via diffusion in SS316 stainless steel under irradiation conditions is crucial. For accurate simulation of this behavior, an interatomic potential encompassing hydrogen, iron, chromium, nickel, molybdenum, carbon, manganese, phosphorus, sulfur, silicon, and nitrogen (H-Fe-Cr-Ni-Mo-C-Mn-P-S-Si-N) would ideally be required. But given the complexity of SS316’s composition and the extremely low concentrations of the elements such as carbon, manganese, phosphorus, sulfur, silicon, and nitrogen, it becomes both impractical and computationally prohibitive to develop a comprehensive potential that includes all these elements. To illustrate, simulating the 0.08 wt% carbon content alone in SS316 would require a simulation cell comprising thousands of atoms, a computationally unfeasible task starting from DFT calculations. However, one carbon atom per 200 atoms in DFT representing a Fe-Cr-Ni-Mo-Mn-C-H system, may be a good starting point to incorporate the effect of carbon inclusion in SS316.

Recent advances have led to the development of a simpler Fe-Cr-H potential¹⁴³ that, while not covering all the complexity of SS316L, can approximate simulations for materials with compositions somewhat similar to SS316L. To better match the specific characteristics of SS316L, it is recommended to extend this potential by incorporating nickel, molybdenum and manganese to create a Fe-Cr-Ni-Mo-Mn-C-H model. This more comprehensive potential would enhance the accuracy of simulations, particularly in studying the effects of irradiation on the tritium transport within the cladding material.

LiAlO₂ pellets

The LiAlO₂ pellets serve as a crucial element within the TPBARs^15,144, where a specified quantity of ⁶Li-enriched material per unit length is incorporated for critical national security to produce tritium at a desired rate. Enrichment of γ-LiAlO₂ with ⁶Li plays a vital role, as this isotope captures neutrons in a light water reactor (LWR) setting, initiating the reaction ({6 atop 3}{Li}+nto {3 atop 1}H+{4 atop 2}{He}+4.8{MeV}). This reaction not only produces tritium but also results in a significant reduction of Li atoms due to irradiation, leading to the formation of a Li-depleted secondary phase, namely LiAl₅O₈^145,146, through the stoichiometric transformation: 5 LiAlO₂ → LiAl₅O₈ + 4 Li (displaced) + 2 O (displaced). This transformation and subsequent phase formation have been corroborated by recent studies¹⁴⁷. The newly formed LiAl₅O₈ phase engages in additional tritium generation that diffuses through the lattice.

With regards to the molecular dynamics simulations for the Li-Al-O system, several interatomic potentials have been developed, including those by Jacob et al.¹⁴⁸, Tsuchihira−Oda−Tanaka (TOT) (2009)¹⁴⁹, and Kuganathan et al.¹⁵⁰. However, Jacob’s potential was deemed unsuitable for our study due to its inability to produce a stable structure in simulations¹⁵¹. Although the TOT potential is well established for the LiAlO₂ system, it is not applicable to both LiAlO₂ and LiAl₅O₈ systems for consistent property comparisons. The fixed ionic charges in the TOT potential (+0.7 for Li, -1.1 for O, and +1.5 for Al) do not ensure charge neutrality in the LiAl₅O₈ system, rendering it non-transferable between these two systems. Roy et al.¹⁵² utilized the potential proposed by Kuganathan et al.¹⁵⁰, which incorporates both short-range interactions (electron-electron repulsion and van der Waals attractions) and long-range Coulombic forces. The ionic charges are allocated as +1 for Li, +3 for Al, and -2 for O (combining core and shell), ensuring charge neutrality across both the LiAlO₂ and LiAl₅O₈ systems. This potential employs the Buckingham model for short-range interactions, characterized by a Pauli repulsion term and a van der Waals attraction term. During cascade simulations where atomic separations may decrease below 0.5 Å, the Buckingham potential is smoothly transitioned to the Ziegler-Biersack-Littmark (ZBL) potential to maintain an effective short-range interaction framework¹³. While most MD works only simulate radiation damage in the LiAlO₂ system^151,153, inclusion of radiation products such as tritium and He have still remained undone. This is due to the lack of reliable interatomic potentials that encompass ({3 atop 1}H) and He. To address this issue, the authors have identified a reactive force field (ReaxFF) developed by Narayanan et al.¹⁵⁴ that encompasses the H-Li-Al-O system. MD simulations by the authors’ group are currently being carried out to simulate tritium behavior in LiAlO₂ using this potential and a good transferability has been noticed. However, for a more accurate analysis, a new ML potential is recommended to be developed for tritium behavior in LiAlO₂ and LiAl₅O₈ during irradiation. The primary requirement is that the potential can handle the high kinetic energy and velocity experienced during irradiation.

Zircaloy-4 getter

The alloy Zr4 consists of Zr (98 wt. %), Sn (1.5%) Fe + Cr (0.37%). However, most modeling efforts using DFT methods have not considered Fe and Cr¹⁵⁵ and have modeled only Sn defective Zr to check the tritium diffusivity in the system. With regards to MD simulations, only HCP Zr has been simulated to assess radiation damage in its pure form without the inclusion of any other metal or even tritium^21,156. While Di et al.¹⁵⁶ focused on the effect of elastic strain on the formation of defects and clusters with a 10 keV PKA to mimic the reactor environment. They used the Mendelev potential¹⁵⁷ which only describes Zr-Zr interactions. Another work by Zhou et al.²¹ investigated up to 80 keV PKAs in hcp-Zr and found that high energy PKAs can create experimental scale vacancy clusters (3 nm size). They used the same Mendelev¹⁵⁷ potential combined with ZBL and also introduced the effect of electronic stopping. Although these works reveal useful insights about Zr, they do not consider the effect of Sn and the behavior of tritium during irradiation. There are certain potentials listed on NIST database that include the H-Zr interactions such as the 2014 potential by Lee et al.¹⁵⁸. However, the potential has not been tested under irradiation environment, and its accuracy cannot be commented upon without sufficient investigation. Nevertheless, this potential can serve as a good starting point to assess the tritium behavior in Zr in the absence of irradiation and at low dose irradiation to start with. Some new machine learning interatomic potentials are under development understand the tritium formation and diffusivity in pure and defective Zircaloy-4 getters¹⁵⁹.

Nickel plating

A Ni plating is deposited around the inner and the outer surface of the Zircaloy-4 getter tube in TPBARs. While several works have investigated radiation damage in pure Ni and its alloys, as listed in Table 1, none of the works have investigated the tritium behavior in Ni with or without irradiation. Tehranchi et al.¹⁶⁰ developed a potential in 2017 to study hydrogen embrittlement in Ni. The same potential may be adopted to study behavior of tritium in Ni with and without irradiation as a starting point. The potential also takes care of grain boundaries which is a major advantage in the case of TPBARs. An alternative potential is also available by Angelo et al.¹⁶¹. Originally developed in 1995, this potential was used to study screw, edge, and Lomer dislocations in Ni and how hydrogen is trapped in the Ni lattice defects. As usual, it’s accuracy for high kinetic energy simulations such as the PKA simulation is yet to be verified. Since the accuracy of the transferability is not tested, ML potential development could prove to be worthwhile for this system.

In summary, this review has examined the landscape of molecular MD simulations to study radiation damage in various materials, including metals, alloys, ceramics, and semiconductors. The intrinsic complexity of radiation damage processes necessitates a detailed understanding of the PKA interactions and the resulting defect cascades. Our analysis has highlighted the significant influence of interatomic potentials on simulation outcomes, emphasizing the criticality of potential selection in accurately predicting defect formation and evolution. Traditional empirical potentials such as the EAM have been widely used due to their availability and computational efficiency. The advent of MLPs has introduced a new paradigm, offering enhanced accuracy by leveraging extensive quantum-mechanical data. These MLPs, including GAP and tabulated GAP (tabGAP), demonstrate superior performance in replicating a wide range of material properties and defect behaviors, albeit with increased computational demands. Our detailed comparison of different materials and potential models underscores the nuanced effects that temperature, PKA energy, and material composition have on radiation damage. For instance, HEAs such as CoCrFeMnNi exhibit higher peak damage, but also enhanced defect recombination compared to pure metals, suggesting a potential route for designing radiation-resistant materials. The inclusion of fission products, such as helium and xenon, in UO₂, further complicates the simulation landscape, requiring specialized potentials to accurately model their interactions and resolution mechanisms. This review also extends to practical applications, providing specific recommendations for interatomic potentials suited for simulating components of TPBARs. For instance, the potential selection for SS 316 L cladding and other materials used in TPBARs is crucial for ensuring the precision of simulations aimed at optimizing tritium production and capture. Recommendations on the basis of potentials for TPBAR components have been provided based on the availability of suitable potentials in the literature. While significant progress has been made in the field of radiation damage simulations, ongoing advancements in machine learning potentials hold promise for even greater accuracy and efficiency. Future work should focus on refining these potentials, expanding their applicability to a broader range of materials and conditions, and integrating multi-scale modeling approaches to bridge the gap between atomic-level simulations and macroscopic material behavior. This holistic approach will be essential for developing next-generation materials capable of withstanding extreme radiation environments, thereby enhancing the safety and performance of nuclear and other high-radiation technologies.