**5. Multiscale approach to modeling radiation damage in materials**

There is a large variety of methods used in condensed matter physics and materials science to study radiation effects in materials, each of them describing a particular aspect of the damage process. **Figure 2** shows a schematic representation of the different time and length scales with the corresponding computational methods that can be applied to study different stages of radiation damage [97, 121, 122]. The very first stage, at the smallest time-length scale, is the electronic stopping regime. For decades, the semiempirical SRIM code discussed in the previous section has been the most widely used tool to calculate electronic stopping power. Nowadays, the electronic stopping power (and the induced electronic excitations in the target) can be described by *ab initio* (parameter-free) methods relying on a realistic description of the electronic and ionic properties of the target system. One of the most accurate methods for treating electronic excitations in materials is timedependent density functional theory (TDDFT) [123] which allows accessing the electronic effects accompanying the ion dynamics. *Ab initio* molecular dynamics (AIMD) [124] based on density functional theory (DFT) [125, 126] can be applied to study point defect formation. At longer time scales and larger length scales, when atomic displacements start to dominate, classical molecular dynamics (MD) and the BCA are usually applied to perform collision cascade simulations. The kinetic Monte Carlo (KMC) [127], the dislocation dynamics (DD) [128], and the finite element method (FEM) [129] are used to study the evolution of defects at the even longer time and larger length scales.

For a complete and accurate description of every aspect of radiation damage, as well as the interplay between them, one has to adopt a combined approach. In recent years, researchers have realized the importance of a multiscale approach to studying radiation damage, as follows from many publications and reviews [121, 122, 130–134]. Each of the methods presented in **Figure 2**, as well the ways of combining them, will be discussed below in the order of increasing complexity. The main focus will be on classical MD, AIMD, and TDDFT, which are fundamental for the description of primary radiation damage at the atomic scale.

### **5.1 Classical molecular dynamics and beyond: the collision cascade**

The most widely used approach in materials science to study the interaction of ions with matter (collision cascades) is MD [135]. MD offers a picture of the ion–ion interaction beyond the linear cascade of the pure BCA by including many-body effects. In MD, atoms are treated as classical particles, and their motion is described by Newtonian dynamics. No electronic effects are thus included.

Cascade simulations need large samples consisting of up to a million atoms (depending on the PKA's energy), which prohibits using parameter-free methods (such as DFT, see Section 5.2) to compute the interatomic forces. Instead, in MD, the forces on atoms are calculated from empirical or semiempirical interatomic potentials (also called force fields) [136–138]. MD with empirical potentials proved to work well for large systems and long time scales [139].

In an MD cascade simulation, the system is usually modeled using periodic boundary conditions, that is, by replicating a small unit cell in all directions. Typically, prior to the cascade simulation itself, a regular MD simulation is done to thermally equilibrate the target system at the desired initial temperature. Then, with the equilibrated configuration, the cascade simulation is initiated by changing the velocity of one of the atoms (the PKA), giving it the desired amount of kinetic energy in the intended direction. The system is then evolved in time as in regular MD, that is, by integrating Newton's equations along with a series of time-steps,

which involves computing the atomic forces, velocities, and positions at each timestep (see Refs. [140, 141] for classical texts on MD). At the end of the cascade simulation, the number of defects is obtained by evaluating the final geometry of the system. Usually, cascade simulations are repeated several times, choosing a different PKA and/or a different direction of the PKA's movement to obtain a statistical average of the number of final defects.

MD has been successfully applied to simulate radiation cascades in a variety of materials [139], from simple metals [142, 143] and compounds [144–146] to complex nanostructures [147], 2D materials [148], and novel multicomponent alloys [149, 150]. MD simulations can afford to access the processes taking place on a relatively long time scale up to ps or even ns which is enough to describe the damage cascade until the thermal spike of the collision has dissipated. Most of the MD codes, however, describe only elastic collisions between atoms and disregard the energy loss mechanisms such as electronic excitation and ionization. The possibility of including electronic excitations is discussed in Section 5.3.

After the primary damage has been formed, defects may continue diffusing, thus annihilating or forming defect clusters. Such processes occur on a much longer time scale, reaching at least seconds, not accessible via regular MD. The problem of simulating a process not accessible in a feasible amount of computational time has motivated the development of several enhanced sampling techniques [151], which in the case of MD simulations of materials have allowed to observe otherwise challenging processes, such as phase transitions.

KMC [127] simulations are commonly used to access long-time effects of radiation in materials [152–155]. KMC is designed to model the time evolution of an atomic system. However, instead of solving the equations of motion, as it is done in MD, the KMC method is based on the assumption that the long-time dynamics of a system consists of diffusive jumps from state to state. Each of the states is treated independently, which makes KMC a very efficient method. The dynamics of the system, that is, the probability of transition from one state to another does not depend on the history of the system. The probability of a state-to-state transition is assigned randomly and the most probable transition is statistically chosen. This allows avoiding the complications related to the choice of interatomic potentials, thus overcoming the time limitations of MD simulations (usually *t*<1 μs) and accessing the macroscopic length scale. KMC is used as an extension of MD to further evolve the damage cascade in time and study the diffusion, accumulation, and annihilation of defects after a collision cascade took place [155].

To further extend the problem into the macro-domain, the DD [128] and FEM [129, 156] methods, based on dividing a geometrical space on a number of finite (non-overlapping) segments, are usually applied. FEM has been used to study the response of a macro-object to external stress in engineering and has also been applied to study the behavior of solids under irradiation by extrapolating the known displacements and evaluating the geometry of a 3D object. DD method allows for calculating the motion of dislocations as well as evaluating the plastic deformation in the material induced by the collective motion of dislocations.

### **5.2** *Ab initio* **molecular dynamics: coupling MD with density functional theory**

AIMD is one of the most important tools in quantum physics and chemistry [157]. In a typical AIMD simulation, it is assumed that the system consists of *N* nuclei and *Ne* electrons for which the Born-Oppenheimer (BO) approximation is applied [158]. The BO approximation implies that the dynamics of the electronic and nuclear subsystems can be treated separately given the fact that the nuclei are much heavier than the electrons and thus the time scales of their motion are very

*Modeling Radiation Damage in Materials Relevant for Exploration and Settlement on the Moon DOI: http://dx.doi.org/10.5772/intechopen.102808*

different. In AIMD, the nuclei are evolved using classical mechanics, while the electronic ground state is adapted to the instantaneous nuclear positions at each step of the dynamics (i.e., the *adiabatic* approximation). The ground-state electronic problem is taken into account through advanced methods, most commonly from DFT [125, 159], a quantum-mechanical method that is used to calculate the electronic structure of many-body systems. Given its importance in describing physical and physicochemical properties of materials, and as it is functional to the understanding of the next section, a brief introduction to the main concepts of DFT is given here.

Practical DFT calculations are based on the Kohn-Sham (KS) formalism [126], which replaces the complex problem of interacting electrons in the standard Schrödinger equation by a problem of non-interacting electrons moving in an effective potential *V*eff:

$$\left\{-\frac{1}{2}\nabla^2 + V\_{\rm eff}([n], \mathbf{r})\right\} \boldsymbol{\upmu}\_i^{\rm KS}(\mathbf{r}) = \varepsilon\_i \boldsymbol{\upmu}\_i^{\rm KS}(\mathbf{r}),\tag{1}$$

where *ε<sup>i</sup>* is the eigenvalues of the KS equations and *ψ*KS *<sup>i</sup>* is the one-electron KS wave functions. Here, the effective potential *V*effð Þ¼ ½ � *n* , **r** *V*extð Þþ **r** *V*Hð Þþ ½ � *n* , **r** *V*xcð Þ ½ � *n* , **r** includes the external potential *V*extð Þ**r** in which the electrons move (i.e., the electron-nuclei Coulomb attraction), the exchange-correlation (XC) potential *V*xcð Þ ½ � *n* , **r** , in which all the many-body effects are included, and the Hartree potential *V*Hð Þ ½ � *n* , **r** which is the electrostatic potential created by the electron density. The solution to self-consistent KS equations Eq. (1), is the exact electron density of the system of interacting electrons, provided that *V*eff is known exactly: *n*ð Þ¼ **r** P*<sup>N</sup> <sup>i</sup>*¼<sup>1</sup> *<sup>ψ</sup>*KS *<sup>i</sup>* ð Þ**<sup>r</sup>** � � � � 2 *:* All the properties of the system (e.g., electronic structure and ground-state energy) can be determined from the electron density, according to the Hohenberg-Kohn theorem [125].

AIMD is used to simulate any physicochemical process where the electronic structure of the system changes significantly or when a detailed description of the structure is needed. A typical example would be the simulation of chemical reactions, where chemical bonds are formed or broken, which cannot be described via classical force fields.

### **5.3 Time-dependent density functional theory for electron dynamics and its coupling to MD**

Although the adiabatic BO approximation is the usual approximation in the methods described above, its applicability is only justified in near-equilibrium situations. However, under ion impact, the electronic subsystem is rapidly driven out of equilibrium.

A realistic description of the dynamics of the electrons in the target during the passage of fast ions can be obtained in the framework of TDDFT which gives access to the electron dynamics out of the electronic ground state. In particular, real-time TDDFT [160] provides a non-perturbative description of the electronic excitations upon an external perturbation and can be combined with the Ehrenfest MD scheme [161], which allows for coupling between electron and ion motion, contrary to the BO picture.

TDDFT consists in solving the time-dependent KS equations [123]:

$$i\hbar\frac{\partial}{\partial t}\psi\_i^{\text{KS}}(\mathbf{r},t) = \left\{-\frac{\hbar^2\nabla^2}{2m} + \hat{V}\_{\text{ext}}(\{R\_I(t)\}) + \hat{V}\_{H\text{XC}}[n(\mathbf{r},t)]\right\}\psi\_i^{\text{KS}}(\mathbf{r},t),\tag{2}$$

where *V*^ *HXC* describes both the electrostatic (Hartree) electron-electron interaction and the quantum-mechanical XC potential, *V*^ *ext* is the potential arising from the ions (both the fast-moving impacting particle and the target atoms), and f g *RI*ð Þ*t* are the atomic positions. The force on the nuclei in Ehrenfest dynamics is defined as

$$F\_I(t) = -\sum\_i \left< \psi\_i^{\text{KS}}(t) | \nabla\_{R\_l} \hat{H}\_\epsilon | \psi\_i^{\text{KS}}(t) \right> \tag{3}$$

where *H*^ *<sup>e</sup>* is the Hamiltonian, that is, the operator on the r.h.s. of Eq. (2). Ehrenfest dynamics is a mean-field method, meaning that the nuclei move on an effective potential energy surface (a mathematical function that describes the energy of the system and whose value depends on the coordinates of all the atoms), which is an average of all adiabatic states involved, weighted by their populations. Other methods exist, in particular, one allowing for electronic transitions, that is, switches between adiabatic states when their population changes [162].

The solution of the time-dependent KS equations in real time can be obtained by applying the so-called time-evolution operator, evolving the KS states in time [123]. The time-step of this propagation must be of the order of attoseconds to describe the fast dynamics of the electrons, in contrast to what occurs in AIMD and MD where the time-step is of the order of femtoseconds. The time-dependent electron density is calculated at each step, from which the total energy of the system is obtained. Knowing the total energy as a function of time, the electronic stopping power can be calculated as *Se* ¼ �*dE=dx*, where *dE* is the energy loss and *dx* is the distance traveled by the projectile inside the target.

Many examples of accurate first-principles calculations of the electronic stopping power are available in the literature [117, 118, 163–168]. Recent studies have demonstrated that electronic excitations (induced by both the primary impacting ion and especially by PKAs and further displaced atoms) affect the cascade evolution [118, 169–171] and thus, they need to be accounted for. The electronic stopping effects can be included in MD cascade simulations through the so-called twotemperature (2T) model [118, 172]. In 2 T-MD, the electrons are included as a thermal bath. Each particle is subject to a friction force representing the electronic stopping and a stochastic force representing the coupling between the vibrational degrees of freedom of the lattice and the electrons. This model considers constant electronic density in the entire system and thus, the electronic stopping power is independent of the crystal direction. Recent studies have extended the 2T model by coupling the electronic and nuclear effects via many-body forces that act in a correlated way. This allowed for the construction of a unified model for ionelectron interactions [170, 171, 173, 174] with a complex energy-exchange process between the ionic and electronic subsystems [174].
