Atomic scale molecular dynamics simulations of radiation damage have been performed on beryllium. Direct threshold displacement simulations along a geodesic projection of directions were used to investigate the directional dependence with a high spatial resolution. It was found that the directionally averaged probability of displacement increases from 0 at 35 eV, with the energy at which there is a 50% chance of a displacement occurring is 70 eV and asymptotically approaching 1 for higher energies. This is, however, strongly directionally dependent with a 50% probability of displacement varying from 35 to 120 eV, with low energy directions corresponding to the nearest neighbour directions. A new kinetic energy dependent expression for the average maximum displacement of an atom as a function of energy is derived which closely matches the simulated data.

## INTRODUCTION

Beryllium (Be) has been chosen as the first wall material for the ITER nuclear fusion reactor,^{1} as its low atomic number and vapour pressure minimise radiative losses from the plasma. Furthermore, Be or high Be compounds have been proposed as neutron multipliers for tritium breeding in the future DEMO reactor.^{2}

Owing to the long history of Be in nuclear applications, its behaviour at elevated temperatures during irradiation has been well characterised experimentally for a thermal neutron flux.^{3–5} More recently, efforts have been undertaken to evaluate its performance under fusion conditions both experimentally^{6–8} and through modelling, whereby the fundamental defect processes, such as the interaction of interstitials, vacancies, impurities, and radiogenic H and He, have been characterized.^{9–12}

One issue that modelling has yet to investigate is the threshold displacement energy (E_{D}) of Be, which is a parameter used in estimations of how many point defects will be created during irradiation. This value may be fed into stochastic Monte-Carlo simulations and approaches based on derivatives of the Kinchin-Pease model, to evaluate the long-term irradiation response of materials.^{13,14}

## METHODOLOGY

Molecular Dynamic (MD) simulations were carried out using the Lammps code.^{15} Through MD, equations of motion for each atom are integrated as a function of time,^{16–18} with forces on each atom calculated from the gradient of the energy, itself calculated from the atom's interaction with its neighbours. This allows the velocities and positions of the atoms to be calculated, with this procedure repeated at discreet time steps, letting the position of the atoms to evolve over time. Canonical (NVT) and isothermal-isobaric (NPT) ensembles are simulated using Berendsen thermostats and barostats.^{19}

Interactions between atoms are described using Embedded Atom Model (EAM)^{20,21} potentials of the form described by Agrawal *et al.*^{22} with a short-range cutoff range of 5 $\xc5$. Note that some constants were misreported in the original paper which has since been corrected with a corrigendum. This potential was chosen as predictions of Be lattice properties corresponding well to experimental data and density functional theory (DFT) simulations, a selection of which are presented in Table I.

Parameter . | Predicted . | Experimental . | DFT . |
---|---|---|---|

E_{c} (eV/atom) | 3.34 | … | 3.32^{a} |

C_{11} (GPa) | 291 | 294^{a} | … |

C_{12} (GPa) | 65.7 | 27^{a} | … |

C_{33} (GPa) | 389.7 | 357^{a} | … |

E_{f}(V) (eV) | 1.26^{b} | … | 0.85^{b} |

E_{f} (Be_{i}) (eV) | 4.69–5.15^{b} | … | 4.20–5.24^{b} |

a ($\xc5$) | 2.372 | 2.286 | 2.27^{b} |

c ($\xc5$) | 3.719 | 3.584 | 3.56^{b} |

c/a | 1.568 | 1.568^{a} | 1.568^{b} |

T_{m} (K) | 1250 | 1550^{c} | … |

Parameter . | Predicted . | Experimental . | DFT . |
---|---|---|---|

E_{c} (eV/atom) | 3.34 | … | 3.32^{a} |

C_{11} (GPa) | 291 | 294^{a} | … |

C_{12} (GPa) | 65.7 | 27^{a} | … |

C_{33} (GPa) | 389.7 | 357^{a} | … |

E_{f}(V) (eV) | 1.26^{b} | … | 0.85^{b} |

E_{f} (Be_{i}) (eV) | 4.69–5.15^{b} | … | 4.20–5.24^{b} |

a ($\xc5$) | 2.372 | 2.286 | 2.27^{b} |

c ($\xc5$) | 3.719 | 3.584 | 3.56^{b} |

c/a | 1.568 | 1.568^{a} | 1.568^{b} |

T_{m} (K) | 1250 | 1550^{c} | … |

The EAM potential takes the form

where E_{i} is the energy of atom i, r_{ij} is the separation between atoms i and j, $\varphi $ is the pair interaction energy, and $F(\rho i)$ is the embedding function, which represents the energy contribution from the interaction of atom i with electron density $\u2009\rho $. The pair potential used is of the Morse form^{25}

where D_{e} is the cohesive energy (0.41246 eV), r_{e} is the nearest neighbour bond length (0.36324 $\xc5$), and α (0.6324 $\xc5$^{−1}) is a parameter related to the curvature about the bottom of the potential well. The density, ρ_{i}, has an exponential form

where A (1.597 e/$\xc5$^{3}) and B (0.49713 $\xc5$^{−1}) are empirically fitted constants. To both the pair potentials and the electron density function, the Voter taper function^{26} is applied to limit the interaction range to a cut-off of 5 $\xc5$ and avoid a discontinuity in the 2nd derivative

The embedding function, $F(\rho i)$, is the Johnson function^{27}

in which F_{0} (2.03930 eV), F_{1} (−12.6178 eV), $\beta $ (0.18752), ρ_{o} (1.0 e/$\xc5$^{3}), and γ (−2.28827) are empirically fitted constants.

Threshold displacement simulations were undertaken by first equilibrating 41 × 41 × 28 Be supercells with dimensions of 94.7 × 94.7 × 100.3 $\xc5$ at a temperature of 300 K for a minimum of 500 ps with timestep 0.001 fs at constant pressure and temperature (NPT) using a Berendsen thermostat and barostat. Periodic boundary conditions were applied. To produce 20 vibrationally different starting configurations, for threshold displacement simulations, equilibration was continued for additional increments of 2 ps to a maximum of 540 ps.

During simulations, primary knock-on energies (PKE) of 5–200 eV in steps of 5 eV are imparted in the desired crystallographic direction to the Primary Knock-on atom (PKA), which is chosen to be the central atom of the cell. Two regions surrounding the PKA, region i and ii, are defined with respect to the application of the thermostat. Region i is a sphere around the PKA with radius optimised to 40 $\xc5$ such that significant displacements of atoms due to the PKA do not extend outside the sphere. To this region, no thermostat is applied. Region ii constitutes the remainder of the supercell and has a Berendsen thermostat applied to remove the excess heat due to the PKA, without directly scaling the atomic velocities in region i.

Following the impact of the PKA, the simulation is run in for 20 000 timesteps of 0.01 fs, 10 000 timesteps of 0.1 fs, and for a further 10 000 timesteps of 1.0 fs, giving a total simulation time of 11.2 ps. This time period was previously shown as sufficient to capture both displacement events and recombination of unstable defects.^{28} It is essential to use such a short timestep during the impact phase due to the high velocity of the PKA and other atoms in the initially damaged region.

The crystallographic directions investigated were those that fall within a 60° clockwise arc from the $[101\xaf0$] direction (as plotted in spherical polar coordinates) and a 90° arc from the (0001) plane to the [0001] direction. Directions within this range were generated based on the vertices of a geodesic projection of the sphere, with a uniform spacing of 6° in both $\phi $ and $\rho $, ensuring results can be directionally averaged to be representative of the overall material similar to Robinson *et al.*^{29} For each direction and energy, 20 replicas were performed. A general schematic of the cell setup is shown in Fig. 1.

## RESULTS AND DISCUSSION

The way the threshold displacement energy, E_{d}, is calculated is by imparting a single atom in a crystal with incrementally increasing PKE in a single direction until it is displaced from its lattice site to another metastable site.^{28} This must be repeated several times for each PKE from slightly different system starting configurations to account for the effects of thermal vibrations on the exact atomic positions at the time of impact. The results of these simulations are presented in Figures 2 and 4.

Figure 2 shows the probability that an atom is displaced at a given energy averaged across all the directions investigated. The lower bound for the threshold displacement energy is PD_{0}, which is the lowest PKE that has a non-zero probability of resulting in a displacement. The model predicts this threshold to be 35 eV for this material. Above 35 eV, the probability of displacement approaches 1 asymptotically, having reached 0.99 for a PKE of 200 eV. The PKE for which there is an average 50% probability that an atom is displaced (PD_{50}) is 70 eV. Both the PD_{0} and PD_{50} have been considered the criteria for threshold displacement energy, with PD_{50} being suggested as more suitable for Monte-Carlo codes such as SRIM.^{29} A further criteria considered are the energy with a threshold displacement probability of 0.1 PD_{0.1}, which is 45 eV. These criteria have been shown to compare favourably with experimental estimations of E_{D} by previous MD studies.^{28} The commonly used Monte-carlo code SRIM uses E_{D} of 25 eV,^{29} while recent molecular dynamics simulations of damage cascades in Be estimate it to be 21 eV.^{30} The large difference is likely attributable to the choice of potential used, as demonstrated for Iron by the work of Nordlund *et al.*^{31} There is also a large difference between the simulated PD_{50} value and that used in SRIM which could have serious implications for the reliability of SRIM calculations in Be. There are no experimental values available to the authors knowledge.

The mean displacement, at the end of the total simulation time (11.2 ps) for a given direction, follows a similar pattern, remaining at approximately 0.45 $\xc5$ (x_{0}) until 35 eV (i.e., PD_{0}). This value is similar to thermal oscillation of the atoms around their lattice sites. Beyond 35 eV, the displacement increases gradually, reaching approximately 4 Å at 200 eV. We can model maximum displacement as a function of PKE by modelling the material as a continuum force field that exerts a drag force on the PKA. We modelled this in two ways: first, with drag as proportional to the momentum of the PKA, and second, as proportional to the kinetic energy of the PKE, giving the two equations describing the drag force on the PKA

where α is the drag coefficient in the momentum model and $\beta $ is the drag coefficient in the kinetic energy dependent model. These can be related to the acceleration through Newton's law to form the differential equations

which given the boundary conditions that at t = 0, $dxdt=2PKEm$, and x = 0, can be solved to give

To calculate the maximum displacement as a function of PKE, at the limit where PKE = PD_{0}, the equations become

These equations are valid when PKE > PD_{0}, below which they have no physical significance. In the absence of a collision event, there will be some thermal motion resulting in a non-zero displacement, x_{0}. Thus, the equations become

These models are compared to the simulated results in Figure 3, with α optimised to 1.06 fs^{−1} and $\beta $ to 1.02. Both models reproduce the simulated results; however, the momentum model (linear drag) overestimates displacement at high PKE, while the energy dependent model is more faithful in this energy range. There is prescience for such an energy dependent correction in the Norgett–Robinson–Torrens (NRT) model which predicts the number of displaced atoms.^{14} The NRT model introduced an energy dependent efficiency term which had the effect of decreasing the number of defects predicted to form with increasing PKA energy.^{14}

The directional dependence of PD_{0} and PD_{50} is examined in Figure 4. By both measures, $ED$ shows strong directional anisotropy, with PD_{0} ranging from 30 to 75 eV and PD_{50} ranging from 35 to 120 eV. Minima in $ED$ are apparent in {1$01\xaf0$} directions, corresponding to the direction of the nearest-neighbour Be in the basal plane. Further minima are observed in {2$1\xaf1\xaf1$}, which correspond to the second nearest Be neighbours out of the plane. Finally, three further minima are apparent surrounding the {0001} directions, which are themselves shallow minima. These relationships hold true for both measures of $ED$. They also support the previous findings of Thomas *et al.*,^{32} who investigated $ED$ in rutile using similar methodology and found that primary knock-on events in the nearest neighbour directions caused a collision sequence resulting in a larger separation of the interstitial—vacancy pair, which was thus more likely to remain stable. Conversely, it was found that glancing-angle collisions resulted in the highest $ED$ as there was little separation between the PKA and vacancy, as well as a larger disordered region which encouraged recombination. Although these results are expected, previous studies have not investigated the directional dependence in such spatial resolution. Such information may be a useful input for binary collision models that take into account the directional dependence of threshold displacement.

## CONCLUSION

We have investigated the threshold displacement energy and its directional dependency using MD in conjunction with classical EAM potentials. It was found that the probability of an atom being permanently displaced increases with primary knock-on energy following a roughly s shaped curve, with a non-zero probability of displacement occurring at 35 eV, and a 50% probability of displacement at 70 eV. The threshold displacement energy was determined to be strongly directionally dependent, with the nearest neighbour directions having a low threshold displacement energy due to the initiation of a collision sequence. Finally, a momentum dependent and an energy dependent model for the maximum displacement of the primary knock-on atom have been developed. The energy dependent model better reproduces the simulated results, especially at the highest displacement energies.

## ACKNOWLEDGMENTS

M.L.J. acknowledges Culham Centre for Fusion Energy for sponsorship. Computing resources were provided by the Imperial College London High Performance Computing Service. Michael Rushton is acknowledged for productive discussions.

## References

_{2}: A molecular dynamics simulation study