# Threshold displacement in the tungsten-carbon system

Matthew Jackson, Paul Fossati, Yan Than Ren, Robin W. Grimes

► To cite this version:

Matthew Jackson, Paul Fossati, Yan Than Ren, Robin W. Grimes. Threshold displacement in the tungsten-carbon system. Philosophical Magazine, 2024, 104 (21), pp.959-976.
10.1080/14786435.2024.2372826.cea-04904972

HAL Id: cea-04904972

https://cea.hal.science/cea-04904972v1

Submitted on 21 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Threshold displacement in the tungsten-carbon system

Matthew Jackson, Paul C. M. Fossati, Yan Ren Than & Robin W. Grimes

**To cite this article:** Matthew Jackson, Paul C. M. Fossati, Yan Ren Than & Robin W. Grimes (2024) Threshold displacement in the tungsten-carbon system, Philosophical Magazine, 104:21, 959-976, DOI: 10.1080/14786435.2024.2372826

**To link to this article:** https://doi.org/10.1080/14786435.2024.2372826

<figure><img src="image_3.png" /></figure>

© 2024 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group

<figure><img src="image_4.png" /></figure>

Published online: 16 Jul 2024.

<figure><img src="image_5.png" /></figure>

Submit your article to this journal

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

Article views: 343

<figure><img src="image_8.png" /></figure>

View related articles

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

View Crossmark data

<figure><img src="image_11.png" /></figure>

CrossMark
OPEN ACCESS

<figure><img src="image_12.png" /></figure>

Check for updates

# Threshold displacement in the tungsten-carbon system

Matthew Jackson<sup>a</sup>, Paul C. M. Fossati<sup>b</sup>, Yan Ren Than<sup>a,c</sup> and Robin W. Grimes<sup>a</sup>

<sup>a</sup>Department of Materials, Imperial College London, South Kensington, London; <sup>b</sup>CEA, Service de Recherche en Corrosion et Comportement des Matériaux, Université Paris-Saclay, Gif-sur-Yvette, France; <sup>c</sup>Singapore Nuclear Research and Safety Institute, Singapore

## ABSTRACT

Threshold displacement energies ($E_d$) of carbon and tungsten in tungsten carbide (WC), $W_2C$, tungsten and diamond are predicted using molecular dynamics. The spatial dependence of $E_d$ is probed by considering a geodesic projection of a symmetrically distinct arc of crystallographic directions for each lattice site. Further, the definition of threshold displacement is explored by making the distinction between atomic displacement ($E_d^{disp}$) and defect formation ($\bar{E}_d^{def}$). Predicted values of $\bar{E}_d^{def}$ compare favourably to experimental observations for tungsten and tungsten carbide. Results confirm that $\bar{E}_d^{disp}$ and $E_d^{disp}$ are strongly structure dependent. Differences between $\bar{E}_d^{disp}$ and $E_d^{def}$ are commensurate with rapid defect recombination within the timeframe of the simulations for some species and structures but not universally. The probability of displacement and defect formation as a function of primary knock-on energy is also reported. Previously developed models for the average displacement of the primary knock-on atom based on kinetic energy and momentum-dependent drag terms are generally found to provide a useful level of approximation. Anisotropy is investigated and results highlight differences due to structures.

## ARTICLE HISTORY

Received 7 October 2023
Accepted 21 June 2024

## KEYWORDS

Radiation damage; tungsten carbide; threshold displacement; molecular dynamics

## 1. Introduction

In the classical description of radiation damage in solids, incident radiation displaces an atom of the material from its lattice site, known as a primary knock-on atom (PKA), which may go on to cause a damage cascade, displacing other atoms through ballistic collisions. The threshold displacement energy, $E_d$, is the minimum energy at which an atom can be displaced from its site. This value can be used to quantify the total damage in a cascade using common

CONTACT Robin W. Grimes

<figure><img src="image_13.png" /></figure>

r.grimes@ic.ac.uk, r.w.grimes@imperial.ac.uk

<figure><img src="image_14.png" /></figure>

Department of Materials,

Imperial College London, South Kensington, SW7 2AZ, London

© 2024 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The terms on which this article has been published allow the posting of the Accepted Manuscript in a repository by the author(s) or with their consent.
<figure><img src="image_15.png" /></figure>

models due to Kinchin-Pease [1], NRT [2] and Greenwood [3]. By extension, if the flux profile of the incident radiation is known then the total number of point defects created in the material can be estimated [2]. This is key in understanding the formation of extended defects, which are responsible for many of the deleterious effects on the macroscopic properties of materials [4,5].

In reality, the probability of displacing an atom with energy E is not a step function from 0 to 1 about E<sub>d</sub>. Thermal oscillations of the atom will add a narrow broadening function relative to typical E<sub>d</sub> for most materials [6,7]. Further, E<sub>d</sub> is strongly directionally dependent (with respect to the crystal structure), with large differences between nearest neighbour directions and low symmetry directions [5–7]. As such, to generate a representative E<sub>d</sub> for a species, the directional dependant probability of displacement, P<sub>d</sub>(θ,φ), must be averaged across all directions:

$$ \bar{P}_d = \frac{\int_0^{2\pi} \int_0^{\pi} P_d(\theta, \varphi) \sin\theta d\theta d\varphi}{\int_0^{2\pi} \int_0^{\pi} \sin\theta d\theta d\varphi} $$

from which $\bar{E}_d$ may be calculated. Two approaches to this are outlined by Robinson et al. [8,9]. In the former, a representative value is calculated by performing displacement simulations across a statistically significant sample of randomly chosen directions. In the latter, an even directional sampling across all angular space is chosen. This has the added advantage that it can accurately probe the directional dependence of E<sub>d</sub>, E<sub>d</sub>(θ,φ).

The comparison of E<sub>d</sub> calculated through molecular dynamics (MD) to experimental results is not simple, and has been explored in detail by Nordlund et al. [5]. Experimentally, E<sub>d</sub> is calculated by irradiating the material with an electron beam and measuring the electrical resistivity, thereby detecting the presence of defects. In atomic-scale simulations, the link between atomic displacement and defect formation is not straightforward. Indeed, atoms can be displaced without creating defects if they swap positions with equivalent neighbours in the crystallographic sense. This leads to the definition of two distinct threshold energies: E<sub>disp</sub>, which is the energy required to displace an atom, and E<sub>def</sub>, which is the energy required to form a defect; which is stable over the time scale of the simulation. In effect, E<sub>disp</sub> is impossible to observe experimentally, and as such experimental values may more likely to be reflected in values of E<sub>def</sub>. However, such behaviour can be observed in MD simulations, given the explicit treatment of atoms. In addition, Nordlund et al. investigated many computational details of displacement simulations [5], which forms a good basis for validation of the models used here.

In this study, the low energy displacement of tungsten and carbon PKAs is investigated and compared across phases in the tungsten-carbon binary system. This system has the advantage of well characterised empirical potentials developed by Juslin and Brenner [10,11]. It also provides a range of crystal structures and chemical environments (see Figure 1) with many applications in the nuclear industry. The structures include diamond, WC, W<sub>2</sub>C, and
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 1** Crystal structures investigated. Dashed atoms denote partial occupancy. Arrows denote [001] direction and length. W C structure from [12].</figcaption></figure>

tungsten. Graphite was also considered, however, initial simulations showed that the carbon potential used in this study, which is successful for WC and W₂C simulations, is unsuitable to simulate graphite due to the lack of a Vander-Waals term.

Tungsten is used as the plasma-facing component of divertors in nuclear fusion test reactors, and may eventually be used as a first wall material [13–15]. Similarly, tungsten carbide is considered for fusion applications [16], particularly as a radiation shield [17,18], and has previously been used as a neutron reflector [4]. As radiation damage may limit the lifetime of materials in these applications, it is vital that the underlying processes are well quantified and characterised. Diamond and W₂C, although unlikely to be used in the nuclear industry, offer a useful basis of comparison for WC and tungsten. Building on the work of Robinson et al., they may be used to assess the extent to which, in these systems, E_d is solely species dependent (as is the default case in SRIM [19]) but depends upon the local chemical or crystallographic environment.

## 2. Methodology

Molecular Dynamic (MD) simulations were carried out using the LAMMPS code [20]. The empirical potential employed is a bond order potential derived by Tersoff et al. for carbon [21], which was later modified by Brenner et al. and then Juslin et al. to describe the entire tungsten-carbon system [10,11]. This potential was selected as it replicates the equilibrium
<figure><img src="image_18.png" /></figure>

properties of the bulk materials investigated in this work well, as shown in Table 1, but further was parametrised by a rigorous fitting regime, which included many non-equilibrium reference states to ensure transferability to different stoichiometries. For example, the work of Juslin et al. [10], reported parameter fitting to alternate (non-equilibrium) structures for tungsten and WC, derived from quantum mechanical calculations based on density functional theory. Defect energies for WC were also computed and compared with previous quantum mechanical simulations, although again, these are static calculations (i.e. effectively at 0 K). The effectiveness of the potential was tested at temperature by predicting the melting point for WC, which was within computational uncertainty of the experimental value, 3049 K. Conversely, defect energies for the tungsten vacancy in tungsten was lower than experimental data. Nevertheless, overall the potential can be reasonably expected to provide a sufficiently good description of the system when far from equilibrium, as would be expected in threshold displacement calculations. Finally, this potential has already been parameterised with a ZBL spline following the procedure laid out by Nordlund et al. [5], ensuring the short-range nuclear repulsion of the two species is adequately described.

The bond order potential of Juslin et al. [10] takes the form:

$$U_{ij} = \sum_{i>j} f_{ij}^{c}(r_{ij}) [V_{ij}^{R}(r_{ij}) - \bar{b}_{ij}V_{ij}^{A}(r_{ij})]$$

where $U_{ij}$ is the total energy of the system, $f_{ij}^{c}(r_{ij})$ a cut-off function, $r_{ij}$ the interatomic separation between atoms i and j, and $\bar{b}_{ij}$ is the bond order term. $V_{ij}^{R}$ and $V_{ij}^{A}$ are repulsive and attractive Morse-like potentials respectively, and take the form:

$$V_{ij}^{R}(r_{ij}) = \frac{D_{0}}{S-1} \exp(-\beta\sqrt{2S}(r - r_{o}))$$

$$V_{ij}^{A}(r_{ij}) = \frac{SD_{0}}{S-1} \exp(-\beta\sqrt{2/S}(r - r_{o}))$$

where $D_0$ is the dissociation energy, $r_o$ the equilibrium bond length, S a fitting parameter and $\beta$ is related to the ground state oscillation frequency.

The cut-off function, $f_{ij}^{c}(r_{ij})$, has the form:

$$f^{c}(r) = \begin{cases} \frac{1}{2} - \frac{1}{2}\sin\left(\frac{\pi(r-R)}{2D}\right), & r \le R - D, \\ 0, & |R - r| \le D, \\ r \ge R + D & \end{cases}$$
<table><caption>**Table 1.** Properties of molecular and bulk compounds in the W-C system predicted using the present bond order potential (BOP).</caption><thead><tr><th colspan="3">W dimer</th><th colspan="3">W-C dimer</th><th colspan="2">C dimer</th></tr><tr><th></th><th>exp.</th><th>DFT</th><th></th><th>exp.</th><th>DFT</th><th></th><th>exp.</th><th>BOP</th></tr></thead><tbody><tr><td>E<sub>c</sub> (eV)</td><td>-2.5</td><td>-2.05</td><td>E<sub>c</sub> (eV)</td><td></td><td>-6.14</td><td>-6.64</td><td>E<sub>c</sub> (eV)</td><td>6.21</td><td>6.00</td></tr><tr><td>r<sub>0</sub> (Å)</td><td>2.2</td><td>1.95</td><td>r<sub>0</sub> (Å)</td><td>1.713</td><td>1.759</td><td>1.905</td><td>r<sub>0</sub> (Å)</td><td>1.243</td><td>1.39</td></tr><tr><td>ω</td><td>337</td><td>-</td><td>ω</td><td>983</td><td>928</td><td>1021</td><td>ω</td><td>1855</td><td>1548</td></tr><tr><td colspan="3">W</td><td colspan="3">Hexagonal WC</td><td colspan="2">C diamond</td></tr><tr><td></td><th>exp.</th><th>DFT</th><th></th><th>exp.</th><th>DFT</th><th></th><th>exp.</th><th>BOP</th></tr><tr><td>E<sub>c</sub> (eV)</td><td>-8.89</td><td>-7.41</td><td>E<sub>c</sub> (eV)</td><td>-16.68</td><td>-15.01</td><td>-16.68</td><td>a (Å)</td><td>3.567</td><td>3.56</td></tr><tr><td>a (Å)</td><td>3.17</td><td>3.16</td><td>a (Å)</td><td>2.907</td><td>2.979</td><td>2.917</td><td>E<sub>c</sub> (eV)</td><td>7.36</td><td>7.36</td></tr><tr><td>B (MPa)</td><td>310</td><td>320</td><td>c/a</td><td>0.97</td><td>0.975</td><td>0.964</td><td>c<sub>11</sub> (MPa)</td><td>1070</td><td>1080</td></tr><tr><td>B' (MPa)</td><td>4.50</td><td>4.20</td><td>B (MPa)</td><td>-</td><td>368</td><td>443</td><td>c<sub>12</sub> (MPa)</td><td>100</td><td>130</td></tr><tr><td>c<sub>11</sub> (MPa)</td><td>522-531</td><td>522</td><td>B' (MPa)</td><td>-</td><td>4.2</td><td>5.1</td><td>c<sub>44</sub> (MPa)</td><td>680</td><td>580</td></tr><tr><td>c<sub>12</sub> (MPa)</td><td>203</td><td>204</td><td>c<sub>11</sub> (MPa)</td><td>720</td><td>651</td><td>710</td><td>E<sub>h</sub>(V) (eV)</td><td>7.2</td><td>7.2</td></tr><tr><td>c<sub>44</sub> (MPa)</td><td>160</td><td>149</td><td>c<sub>33</sub> (MPa)</td><td>972</td><td>887</td><td>896</td><td colspan="2">C graphite</td><td></td></tr><tr><td colspan="3">Rhombohedral W<sub>2</sub>C</td><td>c<sub>12</sub> (MPa)</td><td>254</td><td>183</td><td>224</td><td colspan="2">exp.</td><td>BOP</td></tr><tr><td>E<sub>c</sub> (eV/at)</td><td></td><td>-8.43</td><td>c<sub>13</sub> (MPa)</td><td>267</td><td>189</td><td>305</td><td>a (Å)</td><td>2.464</td><td>2.510</td></tr><tr><td>a (Å)</td><td>3.00</td><td>3.13</td><td>c<sub>44</sub> (MPa)</td><td>328</td><td>-</td><td>267</td><td>c/a</td><td>2.724</td><td>2.714</td></tr><tr><td>c (Å)</td><td>4.77</td><td>4.74</td><td>c<sub>66</sub> (MPa)</td><td>233</td><td>234</td><td>243</td><td>E<sub>c</sub> (eV)</td><td>7.37</td><td>7.380</td></tr></tbody></table>

Note: For the molecules, r₀ is equilibrium bond length, E_c the cohesive energy, and ω the ground state oscillation frequency of a dimer molecule. For the solids, a and c are lattice parameters, B is the bulk modulus, B' is the pressure derivative of the bulk modulus and cᵢⱼ are elastic constants. Experimental and DFT results reproduced from [10] and [11].
<figure><img src="image_19.png" /></figure>

The bond order term, $b_{ij}$ contains three-body and angularity contributions:

$$b_{ij} = (1 + \chi_{ij})^{-1/2},$$

$$\chi_{ij} = \sum_{k(\neq i,j)} f_{ik}^{c}(r_{ik}) g_{ik}(\theta_{ijk}) \omega_{ijk} e^{\alpha_{ijk}(r_{ij}-r_{ik})}$$

where $f_{ik}^c(r_{ik})$ is the cut-off function, $\omega_{ijk}$ and $\alpha_{ijk}$ are fitting parameters, and $g_{ik}(\theta_{ijk})$ is the angular term:

$$g(\theta) = \gamma \left( 1 + \frac{c^2}{d^2} - \frac{c^2}{d^2 + (h + \cos\theta)^2} \right)$$

In which $\gamma$, c, d and h are fitting constants. In addition, at short interatomic separations the repulsive potential $V^R$, is modified with a ZBL term fit to density functional theory data, with the modified potential, $V_{mod}^R$:

$$V_{mod}^R(r) = V^{ZBL}(r)[1 - F(r)] + V^R(r)F(r),$$

$$F(r) = \frac{1}{1 + e^{-b_f(r - r_f)}}$$

where $b_f$ and $r_f$ are species-dependent constants chosen such that the potential at the equilibrium bond length is unmodified, but the potential converges to the ZBL potential at smaller interatomic separations, while ensuring a smooth transition between the two regimes. The empirical constants for W-W, W-C and C-C are reproduced in Table 2.

Supercells for threshold displacement simulations were constructed to have approximately uniform dimensions containing around 5000 atoms, which has previously been shown to be sufficient for such simulations [22]. These were repeated using three-dimensional periodic conditions. This corresponds to an 8 × 8 × 8 supercell containing 4096 atoms for diamond, a 15 × 15 × 15

<table><caption>**Table 2.** Parameter set for the W-C system. C-C parameters from Brenner [11], table reproduced from Juslin et al. [10]. $\alpha_{WWW} = \alpha_{WCW} = 0.45876$, $\omega = 1.0$ for all species.</caption><thead><tr><th>Parameter</th><th>W-W</th><th>W-C</th><th>C-C</th></tr></thead><tbody><tr><td>D<sub>0</sub> (eV)</td><td>5.418 61</td><td>6.64</td><td>6.0</td></tr><tr><td>r<sub>0</sub> (Å)</td><td>2.340 95</td><td>1.905 47</td><td>1.39</td></tr><tr><td>β (Å<sup>-1</sup>)</td><td>1.385 28</td><td>1.803 70</td><td>2.1</td></tr><tr><td>S</td><td>1.927 08</td><td>2.961 49</td><td>1.22</td></tr><tr><td>γ</td><td>1.882 27 x10<sup>-3</sup></td><td>7.285 5 x10<sup>-2</sup></td><td>2.0813x 10<sup>-4</sup></td></tr><tr><td>c</td><td>2.149 69</td><td>1.103 04</td><td>330.0</td></tr><tr><td>d</td><td>0.171 26</td><td>0.330 18</td><td>3.5</td></tr><tr><td>h</td><td>-0.277 80</td><td>0.751 07</td><td>1.0</td></tr><tr><td>R (Å)</td><td>3.50</td><td>2.80</td><td>1.85</td></tr><tr><td>D (Å)</td><td>0.30</td><td>0.20</td><td>0.15</td></tr><tr><td>b<sub>f</sub> (Å)</td><td>12</td><td>7</td><td>8</td></tr><tr><td>r<sub>f</sub> (Å)</td><td>1.3</td><td>1.2</td><td>0.6</td></tr></tbody></table>
supercell containing 6750 atoms for tungsten, and a 14 × 14 × 14 supercell con-
taining 5488 atoms for hexagonal tungsten carbide.

Given the partial occupancy of the carbon sites in the W₂C structure, care
must be taken to ensure a low energy structure is created. To this end, a
10 × 10 × 10 structure containing 2000 W sites was created, with carbon sites
populated according to their probability of being occupied, as observed by
Kurlov et al [12]. This procedure was repeated for 1000 supercells, which
were then relaxed. The 20 lowest energy supercells were then selected to use
as the basis of displacement simulations.

Supercells for all materials were equilibrated at 300 K and zero pressure for a minimum of 50 ps with a timestep of 0.2 fs in NPT conditions, during which temperature and pressure were controlled using the Berendsen thermostat and barostat [23]. Preliminary tests showed that at the small displacement energies we considered, the results were not dependent on the presence of a thermostat at the box boundaries, therefore we did not employ any temperature control scheme after the initial relaxations.

Given the statistical nature of threshold displacement simulations, calcu-
lations must be repeated with different starting configurations. For each
system (except W<sub>2</sub>C, in which 20 configurations have already been generated),
20 starting configurations were generated by continuing the equilibration of the
already equilibrated supercells for additional increments of 1 ps (to a maximum
total equilibration of 70 ps). Threshold displacement simulations were then
performed using the equilibrated cells by impacting a central atom of the
desired species with increments of energy in the specified direction. Energies
between 8 and 160 eV in increments of 4 eV were used for diamond, 20-
200 eV in increments of 5 eV for WC, 4–100 eV in increments of 4 eV for
W<sub>2</sub>C and 20–140 eV in increments of 5 eV for tungsten. Simulations were per-
formed over a total of 12.6 ps using constant volume and energy (NVE) con-
ditions: 20,000 time steps of 0.01 fs during the collisional phase, then 10,000
time step of 0.1 fs and 10,000 time steps of 1 fs during the remainder of the
simulation.

Directions to be investigated were chosen based on a quasi-geodesic projection of directions [24] with spacing in φ and θ of 6°, and covering an angular area including at least double the irreducible symmetry of the structure. This corresponds to the following ranges: 0°-120° and 0°-90° for θ and φ respectively in tungsten carbide and W<sub>2</sub>C, 0°-90° for both θ and φ in tungsten, and 0°-180° and 0°-90° for θ and φ in diamond.

Displacements were detected when any atom had moved and remained at
least 0.6 times the nearest neighbour distance away from its initial position.
This distance is similar to those used previously to study radiation damage-
induced defect formation in beryllium [24] and fluorapatite [25]. The value is
chosen to explicitly exclude displacement due to thermal oscillations [24]. A
defect is said to form if it remains displaced at the end of the simulation
<figure><img src="image_20.png" /></figure>

time. Defects were identified using local environment analysis, as implemented in the coredynamics package [26]. In this method, the tridimensional local average density field (LADF) is used to characterise the local environment of atoms and generate a configuration graph without the need for a reference configuration. Vacancies and interstitials create unique patterns in this configuration graph, which thereby facilitate their identification. Unfortunately, this method was unable to identify defects in the W₂C structure due to the partial occupancy of the carbon sites.

# 3. Results and discussion

## 3.1 Directionally averaged results

Directionally averaged probabilities of defect formation, $\bar{P}_d^{\text{def}}$, and probabilities of displacement, $\bar{P}_d^{\text{disp}}$, as a function of increasing primary knock-on energy, E, are shown in Figure 2. The corresponding threshold defect formation energies, $\bar{E}_d^{\text{def}}$, and displacement energies, $\bar{E}_d^{\text{disp}}$, are calculated by fitting to the model outlined by Robinson et al. [9] and described by the equation.

$$ \bar{P}_d(E) = \begin{cases} 0 & E \le \bar{E}_d \\ \frac{1}{\beta}[E^{\alpha} - \bar{E}_d^{\alpha}] & E \ge \bar{E}_d \end{cases} $$

where $\beta$ and $\alpha$ are fitting parameters. The calculated values of $\bar{E}_d^{\text{def}}$ and $\bar{E}_d^{\text{disp}}$ along with experimental values (where available) are reported in Table 3.

<figure><img src="image_21.png" /><figcaption>**Figure 2.** $\bar{P}_d^{\text{def}}$ (left) and $\bar{P}_d^{\text{disp}}$ (right) as a function of PKA energy in diamond, tungsten, WC and W₂C. Lines are those from the Robinson model [9] fitted to the simulated data points, used to predict the $\bar{E}_d^{\text{def}}$ and $\bar{E}_d^{\text{disp}}$ values reported in Table 3.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Table 3.** Calculated threshold displacement $\bar{E}_d^{disp}$ and threshold defect formation $\bar{E}_d^{def}$ values (eV), with experimental values $\bar{E}_d^{exp}$ or previous molecular dynamic results $\bar{E}_d^{MD}$ where available.</figcaption></figure>

<table><thead><tr><th></th><th>E<sup>def</sup><sub>d</sub></th><th>E<sup>disp</sup><sub>d</sub></th><th>E<sup>exp</sup><sub>d</sub></th><th>E<sup>MD</sup><sub>d</sub></th></tr></thead><tbody><tr><td>Diamond</td><td>21.7 ± 0.4</td><td>12.7 ± 1.0</td><td>37.5-47.6<sup>a</sup>, 35<sup>b</sup></td><td>30<sup>c</sup></td></tr><tr><td>Tungsten</td><td>41.7 ± 0.8</td><td>38.0 ± 0.8</td><td>42<sup>d</sup></td><td>41<sup>e</sup>, 52-68<sup>f</sup></td></tr><tr><td>WC (W PKA)</td><td>47.7 ± 0.9</td><td>41.0 ± 0.9</td><td>42 ± 2<sup>g</sup></td><td>-</td></tr><tr><td>WC (C PKA)</td><td>26.7 ± 1.7</td><td>20.4 ± 1.0</td><td>28 ± 6<sup>g</sup></td><td>-</td></tr><tr><td>W<sub>2</sub>C (W PKA)</td><td>-</td><td>9.3 ± 0.5</td><td>-</td><td>-</td></tr><tr><td>W<sub>2</sub>C (C PKA)</td><td>-</td><td>7.2 ± 0.1</td><td>-</td><td>-</td></tr></tbody></table>

Note: Error is the standard error from multivariant regression of fitting the Robinson model to the spatially averaged data.

a [27] b [28] c [29] d [30] e [31] f [32] g (estimated from $TaC_{0.99}$) [33].

Figure 2 shows that the fitted Robinson model can closely reproduce the simulated $\bar{P}_d^{def}$ and $\bar{P}_d^{disp}$ curves as a function of E. There is, however, some deviation at (low) energies close to where the MD simulations indicate threshold displacement is first observed, with the Robinson model predicting the onset of displacement for C and W in WC higher than MD values.

From Table 3 it is apparent that there is only a modest difference in $\bar{E}_d^{def}$ between a carbon PKA in diamond and WC (21.7 and 26.7 eV respectively). Conversely, there is a significant difference in $\bar{E}_d^{disp}$ for a C PKA in WC compared to a C PKA in diamond and $W_2C$ (20.4 eV compared to 12.7 eV and 7.2 eV respectively). These observations are at variance with the common approximation that $\bar{E}_d$ is solely species dependent, but consistent with being also dependent on the local crystallographic environment. Thus, the current work is consistent with the results of Robinson et al., who found that $\bar{E}_d^{def}$ varied significantly, for both Ti and O, between different phases of $TiO_2$ [8].

Returning to Figure 2, $P_d^{def}$ is similar in diamond and for C PKA's in WC at low PKA energy, E, but at high E, C defect formation becomes significantly more likely in WC. This leads to problems when attempting to use the MD values presented in Figure 2 and in Table 3 in models such as that of Kinchin and Pease (K-P) [1, 7]. The K-P model assumes that the total number of defects generated by a PKA of kinetic energy $E_{PKA}$ is (up to a limiting value and after the second atom is displaced) linearly proportional to the energy required to displace an initial atom, $E_d$ (i.e. the number of displaced atoms, when more than one atom is displaced, is given by $E_{PKA}/E_d$). The values of $\bar{E}_d^{def}$ (reported in Table 3) are the very minimum energy sufficient to cause a defect to be generated (i.e. while defects could be formed, they are not very likely). If these values are used, the KP model rate of C defect production in diamond and WC are almost the same. However, if a value is chosen corresponding to, for example, to a $\bar{P}_d^{def}$ value of say 50% the probability of a C defect being formed (i.e. the rate of success of defect formation) is quite different in diamond and WC. Presumably then, the rate of formation of C defects will be higher in WC than in diamond. An equivalent change in the
<figure><img src="image_23.png" /></figure>

probability of defect formation is predicted when comparing W defect gener-
ation in tungsten and WC.

Some of the differences in $\bar{P}_d^{\text{def}}$ between carbon in different systems may be explained by examining the $\bar{P}_d^{\text{disp}}$ and $\bar{P}_d^{\text{def}}$ curves in Figure 2. The $\bar{P}_d^{\text{disp}}$ curve for a C PKA in WC is consistently modestly greater than $\bar{P}_d^{\text{def}}$ above $\bar{E}_d$, suggesting consistent moderate recombination of defects up to the MD simulation time scale. Comparing these curves for C in diamond, the $\bar{P}_d^{\text{disp}}$ curve raises quickly to saturation whereas the $\bar{P}_d^{\text{def}}$ curve raises much more slowly, consistent with much more effective recombination over this time scale (compared with that in WC). Of course, this comparison between diamond and WC is that of a single-element material and a compound. However, this difference in curves is not predicted for W defects in tungsten and WC, where the curves more closely mirror each other. For a C PKA in $W_2C$, the $\bar{P}_d^{\text{disp}}$ curve is between those of diamond and WC at moderate energies, but it nevertheless gives rise to the lower threshold displacement energy.

$\bar{E}_d^{\text{def}}$ for a W PKA in tungsten and tungsten carbide are 41.7 and 47.7 eV respectively, although interestingly while at low E, $\bar{P}_d^{\text{def}}$ is higher in tungsten, at high E, defect formation is significantly more probable in WC. That $\bar{E}_d^{\text{def}}$ is lower in pure tungsten than in WC may in part be connected to the higher bulk modulus and cohesive energy of WC but of course, also the different crystallography.

Comparing the $\bar{P}_d^{\text{disp}}$ and $\bar{P}_d^{\text{def}}$ curves for tungsten and a W PKA in WC, $\bar{P}_d^{\text{disp}}$ is, of course, always higher than $\bar{P}_d^{\text{def}}$ for both environments, however, the curves are more similar for tungsten carbide, and more different for tungsten. It could be then that recombination is more effective in tungsten, as might be expected due to it being a single element rather than a compound (i.e. in WC, C atoms at C sites may impede W atoms returning to their lattice sites).

The $\bar{P}_d^{\text{disp}}$ curve for a C atom in WC is higher than for a W atom in WC (i.e. for the same energy the probability that a C atom is displaced is greater than a W). Nevertheless, when comparing $\bar{P}_d^{\text{def}}$ curves, for *P* values greater than 0.2, the order is reversed. That is, the $\bar{P}_d^{\text{def}}$ curve for a W PKA in WC is above that for a C PKA (i.e. given the same energy it is more likely a defect will be formed by a W PKA than a C). Again this highlights the necessity to consider carefully which energy to use and at what displacement of defect formation probability in models based around the ideas of Kinchin and Pease.

Comparing the calculated values of $\bar{E}_d$ with experimental and previous MD results ($\bar{E}_d^{\text{MD}}$) in Table 3, $\bar{E}_d^{\text{def}}$ values predicted by the current simulations correspond well for W defects formation in WC and tungsten and for C defects in WC but is lower than that observed experimentally for diamond. However, the $\bar{E}_d^{\text{def}}$ curve for C defect formation in diamond raises significantly more slowly
than other curves, the probability for defect formation is still only ~20% once experimental E<sub>d</sub> values are reached. Nevertheless, another MD study value does predict a higher value, and as such current predictions for diamond should be viewed with more caution.

In addition to examining the probability of displacement and defect formation, another useful measure is the maximum final displacement (x<sub>m</sub>) as a function of E (defined as the displacement of the atom that moved the furthest away from its initial position at the beginning of the simulation). Previously, two models similar to that by Robinson et al. [8] for P&#773;<sup>disp</sup><sub>d</sub> have been developed to predict x<sub>m</sub> based on E<sub>d</sub> for beryllium [24]. One model is based on a momentum-dependent drag term, and the second on a kinetic energy dependent drag term. These models are as follows:

$$ 
\left. 
\begin{array}{l}
x_m = x_0 \\
\\
x_m = \frac{1}{\alpha} \sqrt{\frac{2E}{m}} \left( \sqrt{\frac{\bar{E}_d}{E}} - 1 \right) + x_0 \quad x_m = \frac{2}{\beta} \ln\left(\frac{\bar{E}_d}{E}\right) + x_0
\end{array}
\right\} \begin{array}{l} 0 < E < \bar{E}_d \\ \bar{E}_d < E \end{array}
$$

where α and β are fitting constants, and x₀ is the maximum displacement due solely to thermal oscillations. The simulated x_m as a function of E for each material and PKA species, along with the fitted models are presented in Figure 3.

Both models provide a reasonable description for C PKAs in WC, W₂C and W in tungsten and W₂C, although the momentum-based model tends to be superior, especially in W₂C (although the scatter is large even with this number of simulations). However, for diamond, and to a lesser extent for a tungsten PKA in tungsten carbide, neither model sufficiently describes x_m, as

<figure><img src="image_24.png" /><figcaption>**Figure 3.** Maximum final displacement, x<sub>m</sub>, as a function of PKA energy, E, from MD simulations, with fitted momentum and kinetic energy drag models (dashed and continuous lines respectively).</figcaption></figure>
<figure><img src="image_25.png" /></figure>

there are clear deviations from a simple polynomial or logarithmic function for simulated $x_m$. In both materials, this is manifested by a broad hump around $x_m = 3.2$ Å, suggesting that the isotropic drag model misses some detail at this lower energy range. How important this might be to predictions of damage at higher energies is not clear.

For both W and C PKA's in W₂C, the mean distance travelled rises rapidly with energy beyond $E_d$. This is likely due to the more open W₂C structure, a consequence of the partial occupancy of the C sublattice. As such, a PKA may travel some distance before impacting another atom.

That these drag models can usefully predict the behaviour of these materials and have previously been shown to predict maximum final PKA displacement in beryllium, suggest that they may be adequate for a wide range of materials.

## 3.2 Directionally dependent results

Having examined the directionally averaged threshold displacement results, the directional dependence is now examined. We begin by recalling from section 3.1 that $\bar{E}_d^{\text{def}}$ and $\bar{E}_d^{\text{disp}}$ correspond to the lowest values necessary to form a defect or cause displacement in any direction. In contrast, Figure 4 shows stereographic projections of $E_{d,0}^{\text{def}}$ (still the absolute lowest energy necessary to form a defect or displace a C or W atom) but now as a function of direction. Because it is still a measure of the threshold, so that the probability tends to zero, the subscript d,0 is used. Figure 4 also shows $E_{d,50}^{\text{def}}$ (the threshold, as a function of direction, at which the probability to form a defect is 50%). $E_{d,50}^{\text{disp}}$,

<figure><img src="image_26.png" /><figcaption>**Figure 4.** (001) stereographic projections of tungsten and diamond, showing threshold displacement $\bar{E}_d^{\text{disp}}$ and threshold defect formation $\bar{E}_d^{\text{def}}$ energies as a function of PKA direction, corresponding to the lowest energy $E_{d,0}$ and when the probability of a successful outcome is 50%, $E_{d,50}$ (note different energy scale for $\bar{E}_d^{\text{disp}}$ diamond).</figcaption></figure>
corresponding to the directional dependence for displacement are also shown for tungsten (W) and diamond. It should be noted that the Robinson model is intended to describe $\bar{E}_d^{\text{def}}$ or $\bar{E}_d^{\text{disp}}$ rather than individual directions and cannot be reliably applied to single directions. As such, here $E_d$ is the lowest energy at which displacements (or defect formation) were observed.

In tungsten, both $E_d^{\text{def}}$ and $E_d^{\text{disp}}$ are lowest in and around the <001> family of directions, for both zero and 50% but still moderate in <111>, which are the nearest neighbour directions. The highest $E_d^{\text{def}}$ is in low symmetry directions (e.g. [0.84, 0.20, 0.5]) surrounding the nearest neighbour directions. These high $E_d^{\text{def}}$ directions represent a glancing angle collision with nearest neighbour atoms. In such a collision, the kinetic energy would be distributed depending on the angle between the two atoms, thereby creating two PKA's with a combined energy of the initial PKA, reducing the probability that a permanent defect is formed.

For diamond, starting with $E_d^{\text{disp}}$, Figure 4 indicates that $E_{d,0}^{\text{disp}}$ and $E_{d,50}^{\text{disp}}$, are highest in the two symmetrically equivalent nearest neighbour [111] directions (that are represented in the [100] projection – the other two are in the [100] projection). This is quite different to the prediction for tungsten, where the nearest neighbour displacement values were not unfavourable. As noted from Figure 2, for diamond, values for $E_d^{\text{def}}$ are significantly higher than $E_d^{\text{disp}}$ and in the case of the 50% values show more variation as a function of PKA angle, (note scale difference) whereas for $E_0^{\text{def}}$ there is little angular dependence (again quite different to the prediction for tungsten). What seems to be common between tungsten and diamond for $E_d^{\text{def}}$ is that the least favourable PKA directions are glancing angle collisions.

Figure 5 shows stereographic projections of $E_{d,0}^{\text{def}}$ and $E_{d,50}^{\text{def}}$ in WC for W and C PKAs. Projections of both $E_{d,0}^{\text{def}}$ and $E_{d,50}^{\text{def}}$ for the W PKA show that nearest carbon neighbour directions [211] provide the lowest energy pathway to form a defect, while nearest and second nearest neighbour tungsten directions have higher $E_d^{\text{def}}$. Directions surrounding (glancing) the [001] direction offer the highest $E_d^{\text{def}}$. The angular dependence of $E_{d,0}^{\text{def}}$ and $E_{d,50}^{\text{def}}$ plots are very similar with and $E_{d,50}^{\text{def}}$ values being consistently roughly double.

The case for C PKAs in WC is at first sight somewhat puzzling as there is a marked difference in the high $E_d$ directions for $E_{d,0}^{\text{def}}$ and $E_{d,50}^{\text{def}}$. Directions close to [0001] have the highest $E_{d,0}^{\text{def}}$, but a relatively low $E_{d,50}^{\text{def}}$. In fact, this is due to the scale change and the low $E_{d,50}^{\text{def}}$ value is the same as the high $E_{d,0}^{\text{def}}$ value. What this means is that once the displacement energy is reached and the light carbon atom passes through the tungsten neighbour shell it is very likely to stay displaced and cannot easily return to its initial position. The same is not true for the heavy tungsten atom. This is also reflected in Figure 2, which predicts
<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /><figcaption>**Figure 5.** (0001) stereographic projections of WC, showing for W and C PKA threshold defect formation $\bar{E}_d^{def}$ energies as a function of PKA direction, corresponding to the lowest energy $E_{d,0}$ and when the probability of a successful outcome is 50%, $E_{d,50}$.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Figure 6.** (0001) stereographic projections of W₂C showing for W and C PKAs values of $E_d^{disp}$ for the three carbon sites and single tungsten site. (See figure 2 for site definitions).</figcaption></figure>
that $E_{d,50}^{\text{def}}$ for C is higher than for W but $E_{d,50}^{\text{disp}}$ for C is lower than for W. This is because while for higher energies it is more likely a C atom will be displaced than will W, it is less likely to be permanently displaced to form a defect (stable on this MD time scale). While additional work would be required to be definitive, due to their mass difference, a carbon atom is more likely to bounce back towards its initial position than a tungsten atom after a ballistic collision with one of its neighbours. This is indicative of possible differences between species in compound solids, especially if the constituent elements have very different masses and is worthy of future detailed investigations.

Figure 6 shows stereographic projections of $E_d^{\text{disp}}$ for the carbon sites and tungsten site in $W_2C$. There is significant variation between the three carbon sites, as would be expected given the differences in local coordination. For the C1 site, $E_d^{\text{disp}}$ is lowest in the direction of other C1 sites (which have lowest partial occupancy), moderate in the direction of C2 sites (moderate occupancy) and highest in the nearest neighbour W directions. This is similar for C2 sites, although oddly $E_{d,50}$ is higher in the direction of nearest neighbour carbon sites. For C3, $E_d^{\text{disp}}$ is lowest in the [0001] direction through W trigonal interstices towards the partially occupied C1 site. It is moderate in the nearest neighbour directions of other carbon sites, and highest in the direction of nearest neighbour W sites. W has lowest $E_d^{\text{disp}}$ in directions with nearest neighbour C sites, suggesting that most of the recorded displacements are on the C lattice.

4. Conclusions

The threshold displacement of W or C atoms has been investigated in tungsten, WC and for comparison $W_2C$ and diamond. This is achieved through MD simulations of low energy cascades, employing bond order potentials, initiated by providing a PKA W or C atom with kinetic energy. This may result in defect formation over the time scale of the simulation, with an associated minimum energy $\bar{E}_d^{\text{def}}$, or just displacement of the PKA atom, which in these systems always demands a lower energy $\bar{E}_d^{\text{disp}}$. These energies are averages over a symmetrically distinct arc of crystallographic directions. $\bar{E}_d^{\text{def}}$ values for W in WC and tungsten and for C in WC, are in good agreement with experimental values, but the energy for C in diamond is lower than experiment. Experimental values for $W_2C$ are unknown.

The threshold displacement energy of C is predicted to be different in diamond, WC and $W_2C$, as is the threshold displacement energy of W between tungsten, WC and $W_2C$. This is contrary to the often-used assumption that threshold displacement energy is solely species dependent, but is consistent with the work of Robinson et al. [8] for threshold displacement in $TiO_2$ and
<figure><img src="image_30.png" /></figure>

Tsuchihira et al. in LiAlO₂ [34]. The Robinson model for predicting the prob-
ability of displacement, $\bar{P}_d^{disp}$, as a function of PKA energy was found to be con-
sistent with current results.

Comparing the probability of defect formation $\bar{P}_d^{def}$ as a function of PKA energy identified some interesting behaviour. In particular, for a C PKA in diamond and WC, while $\bar{E}_d^{def}$ values are not dissimilar (21.7 and 26.7 eV) beyond a value of ~0.2, this probability of defect formation increases much more rapidly in WC than in diamond. Also, in WC, for the same PKA energy, the probability that a C atom is displaced is greater than for a W but the probability that this process leads to form a defect is lower for the C PKA than for W. This illustrates that it is not sufficient just to predict the $\bar{E}_d^{disp}$ or even $\bar{E}_d^{def}$ because the likelihood of these processes, as a function of energy, changes differently for different species and in different structures.

For predicting the maximum distance travelled by an atom as a function of its PKA energy, kinetic energy and especially momentum-dependent drag models provide useful descriptions: for a W PKA in tungsten and W₂C, and a C PKA in WC and W₂C. The models are less satisfactory for diamond and for a W PKA in WC, which demand a more complex functional to capture some detail. Nevertheless, given that the drag models were developed and successful for a very different material, beryllium [6], while care must be taken, they are worth assessing for other materials.

Examination of the directional dependence of $\bar{E}_d^{disp}$ and $\bar{E}_d^{def}$ illustrated the strong structural dependence. In general, nearest neighbour directions exhibit high threshold displacement energies, although the highest is typically found at glancing angles to the nearest neighbour.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Funding

Computational facilities and support were provided by High Performance Computing Centre at Imperial College London. Discussions were enabled by the CHASM workshop held at the University of Thessaly.

References

[1] G.H. Kinchin and R.S. Pease, *The displacement of atoms in solids by radiation*. Reports Prog. Phys **18** (1955), pp. 1–51. doi:10.1088/0034-4885/18/1/301

[2] M.J. Norgett, M.T. Robinson, and I.M. Torrens, *A proposed method of calculating displacement dose rates*. Nucl. Eng. Des **33** (1975), pp. 50–54. doi:10.1016/0029-5493(75)90035-7

[3] L.R. Greenwood and R.K. Smither, *SPECTER; neutron damage calculations for materials irradiations*. IJI **3** (1985), pp. 8.
<figure><img src="image_31.png" /></figure>

[4] R.J.M. Konings, T.R. Allen, R.E. Stoller, and S. Yamanaka, *Comprehensive Nuclear Materials*, Elsevier Ltd, Amsterdam, 2012.

[5] K. Nordlund, J. Wallenius, and L. Malerba, *Molecular dynamics simulations of threshold displacement energies in Fe*. Nucl. Instrum. Methods Phys. Res. Sect. B Beam Interact. with Mater. Atoms **246** (2006), pp. 322-332. doi:10.1016/j.nimb.2006.01.003

[6] M.L. Jackson, P.C.M. Fossati, and R.W. Grimes, *Simulations of threshold displacement in beryllium*. J. Appl. Phys **120** (2016), p. 045903. doi:10.1063/1.4958974

[7] Y.R. Than and R.W. Grimes, *Predicting radiation damage in beryllium*. Phil. Mag **101** (2021), pp. 306-325. doi:10.1080/14786435.2020.1834636

[8] M. Robinson, N.A. Marks, and G.R. Lumpkin, *Structural dependence of threshold displacement energies in rutile, anatase and brookite TiO$_2$*. Mater. Chem. Phys **147** (2014), pp. 311-318. doi:10.1016/j.matchemphys.2014.05.006

[9] M. Robinson, N.A. Marks, K.R. Whittle, and G.R. Lumpkin, *Systematic calculation of threshold displacement energies: Case study in rutile*. Phys. Rev. B **85** (2012), pp. 104105. doi:10.1103/PhysRevB.85.104105

[10] N. Juslin, P. Erhart, P. Träskelin, J. Nord, K.O.E. Henriksson, K. Nordlund, E. Salonen, and K. Albe, *Analytical interatomic potential for modeling nonequilibrium processes in the W-C-H system*. J. Appl. Phys **98** (2005), pp. 123520. doi:10.1063/1.2149492

[11] D.W. Brenner, O.A. Shenderova, J.A. Harrison, S.J. Stuart, B. Ni, and S.B. Sinnott, *A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons*. J. Phys. Condens. Matter **14** (2002), pp. 783-802. doi:10.1088/0953-8984/14/4/312

[12] A.S. Kurlov and A.I. Gusev, *Neutron and x-ray diffraction study and symmetry analysis of phase transformations in lower tungsten carbide W 2 C*. Phys. Rev. B **76** (2007), pp. 174115. doi:10.1103/PhysRevB.76.174115

[13] D.M. Duffy, *Fusion power: A challenge for materials science*. Philos. Trans. A. Math. Phys. Eng. Sci **368** (2010), pp. 3315-3328.

[14] Björkas, C. (2009). Interatomic potentials for fusion reactor material simulations.

[15] D. Stork, P. Agostini, J.-L. Boutard, D. Buckthorpe, E. Diegele, S.L. Dudarev, C. English, G. Federici, M.R. Gilbert, S. Gonzalez, A. Ibarra, C. Linsmeier, A.L. Puma, G. Marbach, L.W. Packer, B. Raj, M. Rieth, and M.Q. Tran, *Materials R&D for a timely DEMO: Key findings and recommendations of the EU Roadmap Materials Assessment Group*. Fusion Eng. Des **89** (2014), pp. 1586-1594. doi:10.1016/j.fusengdes.2013.11.007

[16] T. Dash, B.B. Nayak, M. Abhangi, R. Makwana, S. Vala, S. Jakhar, C.V.S. Rao, and T.K. Basu, *Preparation and neutronic studies of tungsten carbide composite*. Fusion Sci. Technol **65** (2014), pp. 241-247. doi:10.13182/FST13-663

[17] C.G. Windsor, J.M. Marshall, J.G. Morgan, J. Fair, G.D.W. Smith, A. Rajczyk-Wryk, and J.M. Tarragó, *Design of cemented tungsten carbide and boride-containing shields for a fusion power plant*. Nucl. Fusion **58** (2018), pp. 076014.

[18] P.A. Burr and S.X. Oliver, *Formation and migration of point defects in tungsten carbide: Unveiling the sluggish bulk self-diffusivity of WC*. J. Eur. Ceram. Soc **39** (2018), 165-172. doi:10.1016/j.jeurceramsoc.2018.10.001.

[19] J.F. Ziegler, J.P. Biersack, and U. Littmark, *The stopping and range of ions in matter*, New York: Pergamon Press, 1985.

[20] S. Plimpton, *Fast parallel algorithms for short-range molecular dynamics*. J. Comput. Phys **117** (1995), pp. 1-19. doi:10.1006/jcph.1995.1039

[21] J. Tersoff, *Empirical interatomic potential for carbon, with applications to amorphous carbon*. Phys. Rev. Lett **61** (1988), pp. 2879-2882. doi:10.1103/PhysRevLett.61.2879
<figure><img src="image_32.png" /></figure>

[22] M.J. Banisalman, S. Park, and T. Oda, *Evaluation of the threshold displacement energy in tungsten by molecular dynamics calculations*. J. Nucl. Mater **495** (2017), pp. 277-284. doi:10.1016/j.jnucmat.2017.08.019

[23] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. DiNola, and J.R. Haak, *Molecular dynamics with coupling to an external bath*. J. Chem. Phys **81** (1984), pp. 3684. doi:10.1063/1.448118

[24] M.L. Jackson, P.C.M. Fossati, and R.W. Grimes, *Simulations of threshold displacement in beryllium*. J. Appl. Phys **120** (2016), pp. 045903. doi:10.1063/1.4958974

[25] E.E. Jay, P.C. Fossati, M.J.D. Rushton, and R.W. Grimes, *Prediction and characterization of radiation damage in fluorapatite*. J. Mater. Chem. A **3** (2015), pp. 1164. doi:10.1039/C4TA01707B

[26] M. Mo, S. Murphy, Z. Chen, P. Fossati, R. Li, Y. Wang, X. Wang, and S. Glenzer, *Visualization of ultrafast melting initiated from radiation-driven defects in solids*. Sci. Adv **5** (2019), pp. 0392.

[27] J. Koike, D.M. Parkin, and T.E. Mitchell, *Displacement threshold energy for type IIa diamond*. Appl. Phys. Lett **60** (1992), pp. 1450-1452. doi:10.1063/1.107267

[28] J.C. Bourgoin and B. Massarani, *Threshold energy for atomic displacements in diamond*. Phys. Rev. B **14** (1976), pp. 3690-3694. doi:10.1103/PhysRevB.14.3690

[29] D. Delgardo and R. Vila, *Statistical molecular-dynamics study of displacement energies in Diamond in diamond*. J. Nucl. Mater **419** (2011), pp. 32-38. doi:10.1016/j.jnucmat.2011.08.035

[30] F. Maury, M. Biget, P. Vajda, A. Lucasson, and P. Lucasson, *Frenkel pair creation and stage I recovery in W crystals irradiated near threshold*. Radiat. Eff **38** (1978), pp. 53-65. doi:10.1080/00337577808233209

[31] C. Björkas, K. Nordlund, and S. Dudarev, *Modelling radiation effects using the ab-initio based tungsten and vanadium potentials*. Nucl. Instrum. Methods Phys. Res. Sect. B Beam Interact. with Mater. Atoms **267** (2009), pp. 3204-3208. doi:10.1016/j.nimb.2009.06.123

[32] Q. Xu, T. Yoshiie, and H.C. Huang, *Molecular dynamics simulation of vacancy diffusion in tungsten induced by irradiation*. Nucl. Instrum. Methods Phys. Res. Sect. B Beam Interact. with Mater. Atoms **206** (2003), pp. 123-126. doi:10.1016/S0168-583X(03)00697-9

[33] D. Gosset, J. Morillo, C. Allison, and C.H. De Novion, *Electron irradiation damage in stoichiometric and substoichiometric tantalum carbides TaCₓ part 1: Threshold displacement energies*. Rad. Effects Def. Solids **118** (1991), pp. 207-224. doi:10.1080/10420159108221360

[34] H. Tsuchihira, T. Oda, and S. Tanaka, *Effects of threshold displacement energy on defect production by displacement cascades in α, β and γ-LiAlO₂*. J. Nucl. Mater **442** (2013), pp. S429-S432. doi:10.1016/j.jnucmat.2013.05.048