Modeling high-energy radiation damage in nuclear and fusion applications

a School of Physics, Queen Mary University of London, Mile End Road, London E1 4NS, UK b School of Physics & Astronomy, University of Southampton Highfield, Southampton SO17 1BJ, UK c Computational Science and Engineering Department, CCLRC Daresbury Laboratory, Daresbury WA44AD, UK d Department of Earth Sciences, University of Cambridge, Downing Street, Cambridge CB2 3EQ, UK e Accelerator Laboratory, University of Helsinki, P.O. Box 43, FIN-00014 Helsinki, Finland


Introduction
Radiation effects are common in nature, with their sources varying from cosmic radiation to decay of isotopes in terrestrial rocks. In addition, a variety of radiation sources are created and used in science and technology. This includes an important area of energy generation in existing nuclear power stations, where kinetic energy of fission products is converted into heat and subsequently electricity. Future fusion reactors, when feasible, will harvest the energy from thermonuclear reactions, mostly from highly intense neutron radiation.
In these and other applications, the energy of emitted particles has a twofold effect. On one hand, this energy is converted into useful energy, by heating the material around the reacting nuclei. On the other hand, this energy also damages this material and degrades its important properties, including mechanical, thermal, transport and so on, to the point that a material might loose its functional purpose. This is currently a central issue for fusion reactors, where the ability of metal structural components to withstand very high neutron fluxes is intensely discussed [1][2][3].
Apart from energy generation in fission and fusion reactors, another important problem is safe encapsulation of highly radioactive nuclear waste. Safe encapsulation of nuclear waste is considered to be one of the pressing issues that modern society faces. It is closely related to public acceptance of nuclear industry and therefore to the future of nuclear power. Regardless of the future of nuclear industry, the amount of accumulated nuclear waste is very large and is steadily growing while no acceptable solution of its safe storage exists.
The highly radioactive nuclear waste mostly comes from spent nuclear fuel, and is radioactive enough to pose a danger to the living environment. The radioactive elements include long-lived actinides. Therefore, the radioactive waste should be safely immobilized over long period of times of up to million years [4][5][6]. The main challenge for safe waste storage is that an encapsulating material should remain an effective immobilization barrier and therefore be resistant to radiation damage from slowly decaying nuclei over this very long time period. Experiments show that radiation damage in some potential waste forms significantly increases diffusion, at which point the material looses its ability to become an effective immobilization barrier. This takes place when local damaged regions connect to form a percolating cluster [7,8].
The problem of resistance to radiation damage is also common in other areas. For example, ion doping can beneficially alter a material's semiconductor properties. On the other hand, the inevitably accompanied structural alteration, if significant, is unwanted, and is therefore referred to as ''radiation damage'' [9,10]. Consequently, semiconductors resistant to radiation damage have important advantages, and are widely studied [11].
In this paper, we discuss radiation damage effects in materials to be used for encapsulation of nuclear waste (waste forms) as well as in fusion reactors. In particular, we focus on modeling these effects using molecular dynamics (MD) simulations. In these simulations, an atom is given high initial velocity, and displaces atoms on its path which, in turn, displace other atoms in the system [10,[12][13][14][15][16][17]. A collection of these atoms is often referred to as ''collision cascade''.
For energies in the range 10-100 keV, the duration of the cascade is typically several ps and its size is several tens of nm. Therefore, MD simulations enable one to directly model the process of creation and evolution of collision cascades on time and length scales not currently accessible experimentally. The insights from these simulations are therefore important for understanding main physical phenomena involved in radiation damage as well as their effect on material properties. As a result, MD simulations have been widely used as a tool complementary to experiments to study and understand radiation damage effects [10,[12][13][14][15][16][17].
Several important limitations exist for MD simulations of radiation damage in nuclear and fusion applications. The main problem is that the energies associated with these processes are very high. Most of the structural damage in waste forms is inflicted by heavy recoils from alpha decay, with associated energies of about 100 keV. The recoil energy is larger in fusion reactors: accelerated by 14 MeV neutrons, iron recoils reach the energy of up to 1 MeV [1,12]. In fission nuclear reactions, the fission product energies are on the order of 50-100 MeV, and fission neutrons can transfer hundreds of keV to atoms in structural components, such as Fe or Zr.
Depending on recoil type, 100 keV cascades require simulating systems with about 10 million atoms whereas 1 MeV cascades require a system on the order of 1 billion atoms in MD simulation cell. Until recently, these systems have remained out of reach for MD simulations. As a result, the energy of simulated collision cascades has typically remained below 100 keV [13][14][15][16][17][18]. At the same time, the need to simulate realistic energy cascades has been particularly emphasized [1]. More generally, the need to simulate length and energy scales that are relevant and appropriate to a particular physical process has been recognized and reiterated [19].
We note that high-energy cascades can branch into smaller subcascades. These can be simulated with lower energy recoils whereas the distribution of sub-cascades can be estimated using approximate schemes such as binary collision approximations [20,21]. However, the energy threshold for cascade branching has not been firmly established because high energy MD simulations have remained out of reach as discussed above. In this sense, high-energy MD simulations of realistic energy recoils are important because they will give a full and detailed picture of what the realistic cascade actually looks like as well as underpin the previous work on cascade branching.
The system size limitation has recently started to lift, primarily owing to two factors: the availability of massive parallel computing facilities with tens of thousands of parallel processors available for a single simulation as well as scalable MD codes. In this paper, we review both the early stages of this process as well as some of the current effort to simulate high energy recoils. We study radiation damage in TiO 2 , a potential waste form as well as other important oxides. We also discuss the ongoing work to simulate high energy radiation damage in Fe, the common material present in the fusion reactor. We visualize the damage and discuss its propagation and evolution, including reversible elastic and irreversible structural damage as well as detailed dynamics of defect atoms. We find that different systems have different ability to recover from radiation damage, and this naturally leads us to the process of resistance to amorphization by radiation damage [22]. We conclude with the current work to extend our simulations to MeV energies and billion-atom systems as well as to account for the electronic energy loss effects using novel algorithms.

Systems and simulation details
We begin with the simulations of radiation damage in rutile TiO 2 , primarily because it is a buffer phase in a multi-phase waste form SYNROC [5,6]. SYNROC consists of several complex titanate oxides designed to chemically accommodate different radioactive species. Therefore, insights from radiation damage in rutile will be useful to understand radiation damage effects in other SYNROC phases.
We start with simulating 100 keV U recoils, corresponding to irradiating conditions in alpha decay. The force field includes short-range potentials as well as Coulomb interactions [23]. We note that the Coulomb interactions slow down the simulation by several times. In their absence, the simulated systems could have been increased in size correspondingly. We subsequently compare the damage in TiO 2 with that in SiO 2 and Al 2 O 3 . For SiO 2 and Al 2 O 3 , we have used empirical potentials from References [24,25], respectively. For Fe, we have used many-body embedded-atom potential [26], modified for better reproduction of several important properties of iron (''M07'' from Ref. [27]). At short distances shorter than 1 Å, interatomic potentials were joined to ZBL short-range repulsive potentials [28].
We have used DL_POLY programme, a general-purpose package designed for large-scale simulations [29,30]. With inherent parallelization based on domain decomposition (DD), it delivers ultimate performance to any high processor count applicable to the problem size with excellent performance of order N and N log N for shortand long-range interactions, respectively. The parallelization design, based on spatial DD, guarantees quasi-equal workload balancing and memory distribution of a problem provided the problem's density distribution is quasi-uniform. The evaluation of the longrange contributions of the Coulomb interactions is based on the Smoothed Particle Mesh Ewald method employing three-dimensional discrete fast Fourier transform operations which benefit from the DD partitioning [29,30]. As a result, DL_POLY demonstrates outstanding scalability for various model systems with the increase of processor count.
Radiation damage simulations were done in a constant energy and volume ensemble at 300 K. This was preceded by equilibration runs in constant pressure ensemble at 300 K. We have implemented several features to handle radiation damage simulations [29,30]. First, we used variable time step to account for faster atomic motion at the beginning of the cascade development and its gradual slowing down at later stages. Second, the MD box boundary layer of thickness of about 10 Å was connected to a constant temperature thermostat at 300 K to emulate the effect of energy dissipation into the sample. Third, we used variable array sizes allocated to parallel processors to account for inhomogeneous density that accompanies the evolution of collision cascades.
MD simulations of radiation damage in oxides were run on the UK HPCx [31] parallel computers, with 256-1024 parallel processors depending on the system. Radiation damage simulations in Fe were run on HECToR National Supercomputing Service [32] using 2048-8192 parallel processors. We have run systems with 2-10 million atoms and 70-100 million atoms in the MD cell for simulations in oxides and iron, respectively. For visualization, we have used Atomeye software [33].

Radiation damage
3.1. Oxides: dynamics of collision cascades, reversible and irreversible deformation and resistance to amorphization Fig. 1 shows the stages of development of the collision cascade due to 100 keV U recoil in TiO 2 , and illustrates one of the main strengths of a MD simulation: detailed visualization of the damage at small lengths and time scales where the cascade develops and relaxes. The frames show approximately 2-3 ps of cascade development. The last frame is at 50 ps, and shows no difference compared to the frames around 3 ps. In Fig. 2 we show a detailed picture of the collision cascade at 50 ps.
An interesting effect is seen in frames 12-13, where a cascade increases in size significantly, followed by substantial size reduction in the final frame. An insight into the origin of this effect comes from the simulation of 40-50 keV recoils [34]. We quantify the amount of radiation damage, and calculate the number of defect atoms. An ''interstitial'' atom is an atom that is further than a certain distance d from any atomic position of the original crystalline lattice. A ''vacancy'' can be defined in a similar way. At each time step, we calculate the number of ''defect'' atoms, N d , the sum of interstitials and vacancies, and show it in Fig. 3. For convenience of comparison of N d in different systems, we choose d = 0.75 Å. Importantly, N d , defined in this way, accounts for both damage production and damage recovery. Damage recovery, or resistance to amorphization by radiation damage, is an important and interesting effect as discussed below in detail.
A notable feature of N d is the large peak of about 10 5 at around 1 ps, followed by the reduction of more than two orders of magnitude at longer times. Interestingly, this peak corresponds to elastic expansion of the lattice around the collision cascade, resulting in a large number of atomic displacements that contribute to the peak of N d . Indeed, we plot all atoms displaced by more than 0.75 Å at 1 ps in Fig. 4. We observe that the highly damaged structure in the core of the cascade is surrounded by the deformed surrounding crystalline lattice, with the pattern of atomic displacements consistent with TiO 2 symmetry.
The origin of this deformation can be discussed as follows. In the absence of radiation damage, atoms vibrate with small amplitudes, and harmonic approximation applies with no volume change. On the other hand, large-scale atomic motions accompanying the development and evolution of the collision cascade considerably widens the distribution of interatomic separations Dr. At large Dr, the harmonic approximation no longer applies, and the potential anharmonicity becomes important. In particular, the decrease of Dr results in a short-range repulsion that is always stronger than the attraction due to the same increase of Dr. Because the expansion carries a smaller energy penalty, large-scale motion of atoms involved in the cascade results in an expansion of the cascade. Consequently, the lattice around the cascade expands, as is seen in Fig. 4. The physical origin of this expansion is therefore related to the anharmonicity of interatomic potential, which also gives rise to the usual expansion in response to temperature increase.
We note that the initial stage of cascade development is often discussed in terms of the ''thermal spike'', invoking the analogy between the liquid motion and atomic motion of the cascade atoms [35]. If the atomic velocities in the cascade are formally converted into temperature, the temperature exceeds the material melting point, suggesting a phenomenological analogy with a molten state. Then, the expansion of the lattice around the  Only atoms that are displaced more than 0.75 Å are shown. This snapshot shows that large number of displaced atoms around the highly damaged core are arranged on the crystalline lattice. This elastic deformation, also seen in frames 12-14 in Fig. 1, is reversible, and corresponds to the large peak of defect atoms in Fig. 3 at 1 ps. Note an interesting ''ray''-shape pattern of atomic displacements that mostly run along the nearest Ti-O neighbor directions, projected on the plane. cascade can be explained because a liquid has, most of the time, lower density compared to a solid. However, the analogy between a liquid and a collision cascade is largely superficial because both structure and dynamics of a liquid essentially differ from those in the collision cascade. On the other hand, our explanation based on the anharmonicity of interatomic potential does not suffer from this problem and offers a microscopic picture of collision cascades instead.
We also note that the damaged structure as well as the lattice around it can expand as a result of the pressure release following the projectile hit that creates high-density regions in the impact area [36]. This expansion is different in origin from the one discussed above: indeed, the proposed mechanism of expansion based on potential anharmonicity does not require density increase, and is based on the widening of interatomic separations only.
Being elastic, the deformation of the lattice around the cascade seen in Fig. 4 is reversible, and disappears after several ps, resulting in the decrease of N d in Fig. 3. The remaining deformation constitutes irreversible structural changes that is ultimately responsible for radiation-induced damage.
In a large class of systems, atomic structure in the damaged regions can be considered as ''X-ray'' amorphous, meaning that extended to system size, the structure shows no Bragg peaks and long-range order [22,34]. Therefore, the ability to recover from radiation-induced structural changes can be called resistance to amorphization by radiation damage.

Predicting resistance to amorphization by radiation damage
Finding the physical origin of resistance to amorphization by radiation damage has a long history, with many models and theories proposed [22]. Basing on extensive analysis of experimental data [22] as well as theory and ab initio simulations [22,37,38], we proposed that the type of interatomic interactions and the nature of chemical bond play a key role in this process. Here, ab initio simulations were used to quantify the electronic density distributions in various systems, and therefore served to predict resistance to amorphization, albeit indirectly.
Importantly, classical MD simulations can directly predict resistance to amorphization by radiation damage. As we have demonstrated for several systems, these simulations reproduce experimental resistance. In Fig. 5, we show collision cascades in SiO 2 and Al 2 O 3 , and observe that the damage is retained in former but recovers in the latter material. This is consistent with experiments which show low and high resistance to amorphization of SiO 2 and Al 2 O 3 , respectively [39].
We note that certain materials (e.g., SiC) that efficiently recover from a single collision cascade can contain a significant concentration of point defects. These materials can be subsequently rendered ''X-ray'' amorphous when defect-rich areas overlap as a result of cascade overlap at high radiation doses.
To quantify this process, we calculate N d for both SiO 2 and Al 2 O 3 , and show it in Fig. 3. We observe that at long times, N d is approximately 10 4 and 10 2 in SiO 2 and Al 2 O 3 , respectively. A detailed inspection of the damaged structure in Al 2 O 3 shows that the damage is mostly due to point defects, in contrast to highly disordered collision cascade in SiO 2 . This is consistent with high and low resistance of SiO 2 and Al 2 O 3 shown in Fig. 5. Fig. 3 shows that in terms of resistance to amorphization, TiO 2 is in between easily amorphizable SiO 2 and highly resistant Al 2 O 3 . This is consistent with experimental behavior: in some experiments, TiO 2 remains crystalline under irradiation at room temperature [40], and in others it can be amorphized, albeit at very high radiation dose [41].
The above results have an important implication: MD simulations can be used to identify future resistant materials as well as study radiation damage in detail. The ability to predict resistance is important, given the current interest in resistant materials for a variety of applications [22], including the safe encapsulation of nuclear waste. Indeed, experiments have shown [7,8] that radiation-induced amorphization can have detrimental effects for the performance of waste forms.
We make two further remarks about predicting resistance to amorphization in MD simulations. First, let us consider the activation barrier, roughly given by the energy required to break the damaged structure enabling damage recovery. It is interesting to ask why empirical potentials, derived at equilibrium, are also able to predict correct activation barriers that govern resistance to amorphization. An interesting insight into this question comes from the observation [34] that potential curvature at equilibrium correlates with the height of the potential energy barrier: flatter (steeper) interaction at equilibrium generally give smaller (larger) activation barriers.
Second, we note that depending on the material and its activation barriers, resistance to amorphization can operate during time scales exceeding those available in MD simulations. However, in many important systems such as zircon, ZrSiO 4 , the size of the collision cascade in MD simulations [7,8,42] agrees with that measured in natural metamict samples that are about 1 billion years old [43]. This implies that damage recovery beyond approximately 10 ps is not significant. In this case, large activation barriers, related to the high energy required to break strong covalent Si-O bonds in order to recover from the inflicted damage imply that the damage is not recoverable in both simulations lasting 10-100 ps and experiments on the time scale of billions of years. Perhaps a surprising result at first glance, it finds an easy explanation as follows.
Let us consider a damaged structure of SiO 2 and calculate how long it takes to re-crystallize at room temperature T r ¼ 300 K. The calculation is similar to the one we used in our recent theory of glass transition [44], where we estimated the time it takes SiO 2 glass to show liquid properties and flow. Let U be the activation energy barrier to re-crystallization related to energy cost to break Si-O bonds. U can be assumed temperature-independent T g 1500 K gives sðT r Þ ¼ 10 67 s, approximately the fourth power of the age of the Universe. The SiO 2 example above illustrates the case of large U and consequently high T g , resulting in very long s and damage recovery times. On the other hand, smaller U can result in the damage recovering on the experimental time scales. Importantly, MD simulations can be reliably used to predict and study resistance to amorphization in the opposite case of small activation barriers and short damage recovery times. This is illustrated in Figs. 3 and 5 showing almost perfect recovery of highly resistant Al 2 O 3 as is seen in the experiments. It is likely that the same fast recovery takes place in other materials such as extremely resistant ZrO 2 , Gd 2 Zr 2 O 7 and so on, proposed to be used in a variety of important applications including immobilization of nuclear waste [22]. This is confirmed by low-energy MD simulations in ZrO 2 [45]. We aim to study these and other materials in our future work. Here, we note that predicting resistant materials on the basis of short-time MD simulations is reliable because no potential long-time recovery effects are involved. Indeed, a material that is resistant in MD simulation will be obviously resistant experimentally.

Iron: current and future simulations
In this section, we briefly discuss our simulations in Fe. This discussion illustrates our current work only, and details will be presented elsewhere.
As discussed in the previous section, we simulated the propagation of U atoms, corresponding to radiation damage in alpha-decay process. In fusion reactors, it is assumed that most of the structural damage originates when a nucleus of the reactor material (typically, Fe) recoils from the scattering of an energetic 14 MeV neutron from the fusion reaction. The energy of Fe recoils can reach 1 MeV [1,12], significantly higher than the energy of recoils in alpha-decay. This increases the required of the MD simulation box substantially. Further, cross-section of Fe atom is significantly smaller than that of U atom, reflected by the smaller effective atomic radius simulated by ZBL repulsive potential. Consequently, Fe recoil looses energy less efficiently in terms of the number of displaced atoms, pushing the required MD cell to even larger sizes. Finally, reliable potentials for Fe that we use contain many-body terms, which considerably slows down the simulations. The combination of these factors makes highenergy radiation damage simulations in Fe a challenge even for currently existing massive parallel computing facilities.
In Fig. 6, we show the initial stage of radiation damage created by 200 keV Fe recoil in iron. In addition to damage in fusion reactors, this simulation also represents radiation damage from the iron recoil from a fast neutron scattering event in a fission reactor. The system size is about 100 million atoms, and the box size is approximately 1000 Å.
We are currently extending this work in two important respects. First, we are in the process of including the electronic energy loss in the integration algorithm [46,47]. The losses of energy related to electronic excitations become more pronounced as recoil energy increases, and become essential at high MeV energies. They should therefore be particularly important for collision cascades in fusion reactors and other important applications.
Second, by increasing the number of parallel HECToR processors [32] to 16,384 or perhaps beyond, we aim to simulate MD systems of up to 1 billion atoms. This will enable us to simulate collision cascades with energies close to MeV or above, depending on the type of material and recoil. Apart from accessing the energy and length scales relevant for several important processes, including in fusion as well as fission reactions, this may bring new effects and physics related to radiation damage.

Conclusions
In summary, we discussed our recent and current work aimed at simulating high energy radiation damage in materials relevant for encapsulation of nuclear waste as well as materials to be used in advanced fission and fusion reactors, including several important oxides and iron. We studied various stages of evolution and relaxation of high-energy collision cascades, and identified reversible elastic and irreversible inelastic structural changes. We explained the elastic expansion of the lattice around the cascade in terms of anharmonicity of interatomic interactions. The remaining irreversible structural change is related to resistance to amorphization by radiation damage. This resistance was quantified by the number of remaining defect atoms in the damaged structure. We discussed how MD simulations can predict experimental resistance to amorphization, including the important case of highly resistant materials. Finally, we discussed our current work to simulate radiation damage of MeV energies and system sizes of the order of billion atoms using massive parallel computing facilities.