# Atomistic Study of Irradiation-Induced Plastic and Lattice Strain in
Tungsten

Jintong Wu  Department of Physics, University of Helsinki, Post-office box 43,
FIN-00014 University of Helsinki, Finland  Daniel R. Mason  UK Atomic Energy
Authority, Culham Science Centre, Oxfordshire OX14 3DB, UK  Fredric Granberg
<EMAIL_ADDRESS>Department of Physics, University of Helsinki,
Post-office box 43, FIN-00014 University of Helsinki, Finland

###### Abstract

We demonstrate a practical way to perform decomposition of the elasto-plastic
deformation directly from atomistic simulation snapshots. Through molecular
dynamics simulations on a large single crystal, we elucidate the intricate
process of converting plastic strain, atomic strain, and rigid rotation during
irradiation. Our study highlights how prismatic dislocation loops act as
initiators of plastic strain effects in heavily irradiated metals, resulting
in experimentally measurable alterations in lattice strain. We show the onset
of plastic strain starts to emerge at high dose, leading to the spontaneous
emergence of dislocation creep and irradiation-induced lattice swelling. This
phenomenon arises from the agglomeration of dislocation loops into a
dislocation network. Furthermore, our numerical framework enables us to
categorize the plastic transformation into two distinct types: pure slip
events and slip events accompanied by lattice swelling. The latter type is
particularly responsible for the observed divergence in interstitial and
vacancy counts, and also impacts the behavior of dislocations, potentially
activating non-conventional slip systems.

Tungsten; Irradiation; Molecular dynamics; Strain evolution; Lattice swelling

##  I Introduction

The consequences of irradiation in materials, mainly seen as microstructural
evolution, originates from the interaction between the incident particles and
the lattice atoms of the material. Among them, radiation damage produced
during the collision cascade process leads to swelling, embrittlement, and
creep. These mechanical properties determine the lifetime of nuclear plant
parts. Tungsten (W) has been chosen to play a pivotal role in fusion test
reactors under construction, such as the International Thermonuclear
Experimental Reactor (ITER) [1]. Here it will be used as the material facing
the plasma inside the reactor, owing to its high atomic mass, leading to a low
sputtering yield [2], high melting temperature and high thermal conductivity
[3, 4, 5].

Irradiation induces significant changes in both the microstructure and
mechanical properties of materials [6, 7, 8, 9]. The response of tungsten to
irradiation has been extensively studied both experimentally and
computationally. On the experimental side, Reza et al. [10] used transient
grating spectroscopy (TGS) and found that 20 MeV self-ion irradiation showed
defect saturation at doses between 0.06 and 0.1 dpa. Another positron
annihilation spectroscopy study by Hollingsworth et al. [11] showed that
vacancies appear as small clusters at low doses at room temperature. Recent
research suggests that the accumulation of small vacancy clusters due to
irradiation can eventually lead to the formation of macroscopic voids,
resulting in observable swelling of the material [12]. Another study
elucidates that this swelling phenomenon primarily arises from the
agglomeration of self-interstitial defects into dislocation loops and extended
dislocation networks. These dislocations are produced during collision
cascades initiated by, for example, neutron irradiation, leading to the
accumulation of uncompensated relaxation volumes [13]. Radiation promotes the
glide of dislocations at high doses, thereby further extending the dislocation
network. This process involves the thermally activated motion of screw
dislocations, particularly evident in BCC metals [14]. Understanding the
underlying mechanism of plastic flow at high-dose irradiation is crucial to
gaining insights into the macroscopic swelling mechanism. On the simulation
side, high-dose simulations are computationally heavy and more efficient
methods have been proposed to generate characteristic high-dose
microstructures. Derlet et al. [15] proposed the creation relaxation algorithm
(CRA) to quickly reach high dose levels without full MD cascades. Defect
densities obtained by the CRA are found to be higher than experimental values
due to the lack of thermal effects [16]. To address this, combination of CRA
and cascade simulation, known as cascade annealing (CA), have been used and
showed good agreement with full MD simulations at similar doses [17, 18]. In
the high-dose irradiation simulation scenario, particularly in the cell
following CRA for sped-up dose accumulation, the atomic structure becomes
highly disordered, and the complex defect network hinders our ability to
accurately characterize the deformation process within the atomic system
during the CA process.

Nevertheless, understanding the intricate and nonlinear dynamics of the damage
microstructure at high irradiation doses is a multifaceted task, encompassing
various length and timescales, while also being influenced by exposure and
environmental conditions [19]. Consequently, creating a concise atomistic
model that accurately encompasses all relevant factors at the appropriate
scales remains a formidable challenge. The detectable fluctuations in strains
and stresses within irradiated materials provide an avenue for directly
confirming the accuracy of real-space simulations. Elasticity equations serve
as the cornerstone for connecting atomic-scale defects to macroscopic strains.
In classical MD simulations, elasto-plastic deformation is inherently coupled,
and quantifying the magnitude of elasto-plastic deformation to gain a deeper
understanding of the underlying defect evolution mechanism remains
challenging. Earlier works [20, 21, 22, 23, 24, 25, 26, 27, 28, 29] explored
the plastic deformation of nanoscale metal crystals at low temperatures using
MD simulations, but these qualitative studies lacked quantitative magnitude of
plastic deformation. Vo et al. [30] proposed a decomposition algorithm
primarily based on dislocation motion-induced plastic deformation. However, it
is only applicable to one-dimensional configurations. In cases of complex
dislocation motion, such as high defect density cases, the interaction between
dislocations often renders the decomposition results inaccurate. Another
previous study provided a decomposition method [31], based on the virtual
intermediate configuration concept. However, if the identification algorithm
fails to accurately determine the local atomic structure, it cannot perform
elasto-plastic decomposition. For instance, a cell exposed to high-dose
irradiation can cause the algorithm to fail.

In addition to advances in simulation methodology, there have been recent
advances in microstructural characterization. The Wigner-Seitz (W-S) method
[32] has been used extensively to identify and characterize isolated defects
embedded in a perfect reference crystal, but fails when the reference lattice
itself can evolve, or when large local displacements comparable with the
lattice spacing exist [33]. Recently Mason et al. [17] developed a method
based on detecting the Wigner-Seitz cell _locally_ , and using this to detect
void isosurfaces. We refer to this method as Void Isosurface Detection (VID).
Machine Learning analyses have also been developed to detect defect features
related to vacancies and interstitials [34, 35].

In summary, achieving a clear, unambiguous, and consistent interpretation of
data from ion-irradiated materials, particularly in the high-dose regime,
remains a significant challenge. Existing experimental interpretation models
rely on kinetic equations that involve numerous parameters but fail to account
for the microscopic, fluctuating stresses and strains driving defect
interactions at the nanoscale, as evidenced by prior studies [36, 37, 38].
Therefore, leveraging recent algorithmic advancements, we attempted to answer
the following key research inquiries in this paper:

  1. 1.

How can we construct a quantitative explanatory model for irradiation effects
at high doses? Can we effectively employ the microscopic fluctuating strain
model as a bridge connecting atomic-scale defects to macroscopic strains [13]?
Additionally, can various defect detection methods accurately signify the
accumulation or abrupt release of internal loading resulting from high-dose
irradiation at the atomic level?

  2. 2.

In simulations, lattice strain measurements can be derived from atomic
positions via diffraction pattern peak positions. However, performing this
calculation for overlapping cascade MD simulations involving tens of thousands
of frames is prohibitively resource-intensive. Are there alternative
approaches to acquire plastic and lattice strain information, directly from
real-space simulations, even in the high dose regimes?

  3. 3.

Under different dosage conditions, how do dislocation loops behave, and do
they induce plastic deformation? If so, does the slip only occurs on certain
planes? Furthermore, can we further categorize the specific type of plastic
transformation involved?

In this article we provide a quantitative interpretation of the role played by
prismatic dislocation loops as initiators of plastic strain effects in heavily
irradiated metals. This phenomenon leads to experimentally measurable
alterations in lattice strain. First, by employing different defect
characterization techniques, we introduce the “anomalies” observed in high-
dose conditions (as detailed in section III.1). To find the mechanisms behind
these phenomena at high doses and provide answers to pertinent questions, we
show how atomic and plastic strain can be identified from snapshots of
molecular dynamics simulations (section II.2). Importantly, we do this without
tracking the relative positions of the atoms, as irradiation mixing makes this
difficult. In section III.2.1 we show how the accumulation of lattice defects
at low dose leads to an increasing homogeneous lattice strain, and how
individual defects can reorient themselves spontaneously to reduce elastic
energy. Finally in section III.2.2 we show how plastic strain starts to emerge
at high dose, and so dislocation creep and irradiation-induced swelling both
arise spontaneously through the agglomeration of dislocation loops into a
dislocation network. We conclude that simulated irradiation shows both lattice
and plastic strain responding to the evolution of defects on molecular
dynamics timescales. This demonstrates unambiguously that a high density of
irradiation induced defects self-organises almost instantaneously to reduce
elastic energy density by dislocation-induced plastic slip.

##  II Simulation Methods

###  II.1 Irradiation simulations

All simulations were conducted with LAMMPS [39], implemented with an adaptive
timestep [40, 32, 41]. The dimensions of the cell were $120\times 120\times
120$ conventional BCC unit cells, resulting in 3.456 million atoms. Periodic
boundary conditions were applied in all dimensions. A friction force was
applied to atoms according to the energy loss table provided by “Stopping and
Range of Ions in Matter” (SRIM) [42] code to include effects of electronic
stopping power, when the kinetic energy of an atom exceeds 10 eV [43]. We used
the embedded atom method (EAM) potential by Ackland and Thetford [44] with the
short-range modified by Zhong et al. [45], AT-ZN. This potential has
previously shown a good agreement with experimental results [46].

For full MD cascade simulations, a 30 keV PKA was introduced at the centre of
the system, in a random direction. Most of the cell was able to evolve without
any constraints, except a two nm thick layer at the borders which were
thermally controlled to 300 K. No pressure control was active during this
simulation. Each cascade was followed for 30 ps, followed by a relaxation
period of 10 ps with temperature control at 300 K on all atoms, and stress-
strain boundary conditions
$\sigma_{x}=\sigma_{y}=\sigma_{z}=\varepsilon_{xy}=\varepsilon_{yz}=\varepsilon_{zx}=0$.
Thereafter, the cell was shifted before the next cascade event, in order to
obtain a homogeneous irradiation. This procedure was repeated 4000 times to
achieve the dose of $\sim$0.11 dpa (with a threshold displacement energy of 90
eV) according to the classical Norgett-Robinson-Torrens displacements per atom
(NRT-dpa) model [47, 48, 49]. Following Ref. [17], we find the canonical DPA
level, cdpa, defined as the expected atomic fraction of vacancies produced.
This value was obtained from the first 30 overlapping cascades to be
$\mathrm{cdpa}=4.2\times 10^{-6}N_{\mathrm{casc}}$.

The CRA method [15] directly inserts FPs to quickly construct a damaged
system. For this work, we randomly inserted 1000 interstitials and deleted
1000 atoms per step, and relaxed the system with the conjugate gradient method
(CG), similarly to previous studies [15, 50]. This was repeated enough times
to achieve doses of 0.01, 0.0143, 0.03, 0.1 and 0.2 cdpa. As the CRA method
does not include either thermal annealing or cascade annealing (CA), to
produce a high dose microstructure comparable with full MD simulations we
subsequently performed CA with 1600 30 keV cascades (0.0078 cdpa) on the CRA
simulation boxes. This combination of CRA+MD has been shown to produce
simulated irradiation microstructures almost indistinguishable from MD alone
[17, 18]. A long thermal annealing was also investigated, however, this method
did not produce comparable results to full MD, like the CA did. More details
about the results and confirmation of convergence in the combination of CRA+MD
used in this work can be found in the Supplementary Information (SI).

All quantitative simulation results presented are the average of three
different simulation runs carried out.

###  II.2 Atomic vs plastic strain in a single crystal

In this section we derive the relationship between lattice strain and plastic
strain in a single crystal simulation. The _total_ strain,
$\underline{\underline{\mathbf{\varepsilon}}}$ , is well-defined, being the
change in dimensions of the supercell. But this strain can be decomposed into
different contributing factors. The _elastic_ strain,
$\underline{\underline{\mathbf{e}}}$, is related by Hooke’s law to the stress
on the boundary via the elastic compliance tensor, ie
$\underline{\underline{\mathbf{e}}}=\underline{\underline{\mathbf{C}}}^{-1}\underline{\underline{\mathbf{\sigma}}}$.
In Mura’s formulation [51], the total strain
$\underline{\underline{\mathbf{\varepsilon}}}=\underline{\underline{\mathbf{e}}}+\underline{\underline{\mathbf{\varepsilon}}}^{\star}$,
where the second term $\underline{\underline{\mathbf{\varepsilon}}}^{\star}$,
known as the _eigenstrain_ , is the _stress-free_ strain in the body due to
defects. Here we are treating the atoms explicitly, and so can make a further
distinction. A part of the strain is visible to the atoms, visible in the
location of the x-ray diffraction peaks, comes from the eigenstrains arising
from the irradiation induced defects. This strain affects the evolution of
small irradiation-induced defects, as they have some freedom to rotate to
compensate the elastic energy density. But a second part of the strain comes
from plastic events, such as dislocation slip or creep, which after having
occurred leave the crystal lattice essentially unaffected. These plastic
strains appear in the diffraction pattern in the detailed fine-structure of
each individual diffraction peak, not in its location.

Consider a periodic supercell with repeat vectors $\mathbf{A}_{1}$,
$\mathbf{A}_{2}$, $\mathbf{A}_{3}$, containing a single crystal with primitive
lattice vectors $\mathbf{b}_{1}$, $\mathbf{b}_{2}$, $\mathbf{b}_{3}$. We can
compactly write these vectors as the matrix
$\underline{\underline{\mathbf{A}}}$, whose columns are $\mathbf{A}_{i}$, and
the matrix $\underline{\underline{\mathbf{B}}}$ whose columns are
$\mathbf{b}_{i}$. Then for the lattice to be continuous across the supercell
periodic boundary requires [17]

$\underline{\underline{\mathbf{A}}}=\underline{\underline{\mathbf{B}}}\,\underline{\underline{\mathbf{N}}},$
(1)

where $\underline{\underline{\mathbf{N}}}$ is a $3\times 3$ matrix of integers
representing the number of lattice repeats. The number of lattice sites in the
periodic supercell is
$N_{\mathrm{latt}}=\mathrm{det}[\underline{\underline{\mathbf{N}}}]\times
N_{\mathrm{motif}}$, where $N_{\mathrm{motif}}$ is the number of atoms in the
motif of the primitive lattice.

The two matrices $\underline{\underline{\mathbf{B}}}$ and
$\underline{\underline{\mathbf{N}}}$ define a _reference lattice_ for the
simulation, compatible with the simulation box size and shape.

Now consider the box of atoms after some simulated irradiation. If we are
running with a general set of stress- and strain-boundary conditions, the
periodic supercell now has repeat vectors
$\underline{\underline{\mathbf{A}}}^{\prime}$. The change in the box can be
written
$\underline{\underline{\mathbf{A}}}^{\prime}=\underline{\underline{\mathbf{F}}}\,\underline{\underline{\mathbf{A}}}$,
where the deformation is given by the product of linear atomic and plastic
deformations, denoted with subscripts a and p, respectively, i.e.
$\underline{\underline{\mathbf{F}}}=\underline{\underline{\mathbf{F}}}_{a}\,\underline{\underline{\mathbf{F}}}_{p}$.
The primitive cell can be atomically strained, and can also be rotated 111To
see that the primitive cell can be rotated, we use proof by _reductio ad
absurdum_ : The most extreme irradiation disordering could be taking out all
the atoms, replacing them in amorphous positions, and annealing back to a
crystal lattice. For no rotation to be true it is necessary that the exact
same lattice orientation is always recovered from amorphisation + annealing.,
so the lattice vectors can also change. We assume here that no second phase
emerges, and a single grain is present after irradiation. If we can compute
the best fit $\underline{\underline{\mathbf{B}}}^{\prime}$ to the atomic
positions after irradiation, then we can write
$\underline{\underline{\mathbf{B}}}^{\prime}=\underline{\underline{\mathbf{F}}}_{a}\,\underline{\underline{\mathbf{R}}}\,\underline{\underline{\mathbf{B}}}$,
where $\underline{\underline{\mathbf{R}}}$ is a rotation matrix. We can use
the method of polar decomposition to find the atomic strain

$\displaystyle\underline{\underline{\mathbf{F}}}_{a}\,\underline{\underline{\mathbf{F}}}_{a}^{T}$
$\displaystyle=$
$\displaystyle\left(1+\underline{\underline{\mathbf{\varepsilon}}}_{a}\right)\left(1+\underline{\underline{\mathbf{\varepsilon}}}_{a}\right)^{T}$
$\displaystyle=$
$\displaystyle\left(\underline{\underline{\mathbf{B}}}^{\prime}\underline{\underline{\mathbf{B}}}^{-1}\underline{\underline{\mathbf{R}}}^{T}\right)\left(\underline{\underline{\mathbf{B}}}^{\prime}\underline{\underline{\mathbf{B}}}^{-1}\underline{\underline{\mathbf{R}}}^{T}\right)^{T},$
$\displaystyle\underline{\underline{\mathbf{\varepsilon}}}_{a}$
$\displaystyle\approx$
$\displaystyle\frac{1}{2}\left(\left(\underline{\underline{\mathbf{B}}}^{\prime}\underline{\underline{\mathbf{B}}}^{-1}\right)\left(\underline{\underline{\mathbf{B}}}^{\prime}\underline{\underline{\mathbf{B}}}^{-1}\right)^{T}-\underline{\underline{\mathbf{1}}}\right).$
(2)

where to get the last line we have taken a linear approximation. From the
arguments above, the matrix of lattice repeats after irradiation must solve
$\underline{\underline{\mathbf{A}}}^{\prime}=\underline{\underline{\mathbf{B}}}^{\prime}\underline{\underline{\mathbf{N}}}^{\prime}$.
As $\underline{\underline{\mathbf{A}}}^{\prime}$ is read from the size of the
MD simulation box, we can rearrange to give the plastic strain

$\underline{\underline{\mathbf{F}}}_{p}=\underline{\underline{\mathbf{R}}}\,\underline{\underline{\mathbf{B}}}\,\underline{\underline{\mathbf{B}}}^{\prime-1}\underline{\underline{\mathbf{A}}}^{\prime}\,\underline{\underline{\mathbf{A}}}^{-1},$
(3)

or as a function of the lattice repeats,

$\underline{\underline{\mathbf{F}}}_{p}=\underline{\underline{\mathbf{R}}}\,\underline{\underline{\mathbf{B}}}\,\underline{\underline{\mathbf{N}}}^{\prime}\,\underline{\underline{\mathbf{N}}}^{-1}\underline{\underline{\mathbf{B}}}^{-1}.$
(4)

Using polar decomposition we can find the plastic strain. As that the matrix
product
$\underline{\underline{\mathbf{N}}}^{\prime}\,\underline{\underline{\mathbf{N}}}^{-1}$
in Eq. 4 is not symmetric, we write
$\underline{\underline{\mathbf{F}}}_{p}=\underline{\underline{\mathbf{R}}}^{\prime}\left(1+\underline{\underline{\mathbf{\varepsilon}}}_{p}\right)$,
where $\underline{\underline{\mathbf{\varepsilon}}}_{p}$ is the (symmetric)
plastic strain tensor. Then

$\displaystyle\underline{\underline{\mathbf{F}}}_{p}^{T}\,\underline{\underline{\mathbf{F}}}_{p}$
$\displaystyle=$
$\displaystyle\left(1+\underline{\underline{\mathbf{\varepsilon}}}_{p}\right)^{T}\left(1+\underline{\underline{\mathbf{\varepsilon}}}_{p}\right)$
$\displaystyle=$
$\displaystyle\left(\underline{\underline{\mathbf{R}}}\,\underline{\underline{\mathbf{B}}}\,\underline{\underline{\mathbf{B}}}^{\prime-1}\underline{\underline{\mathbf{A}}}^{\prime}\,\underline{\underline{\mathbf{A}}}^{-1}\right)^{T}\left(\underline{\underline{\mathbf{R}}}\,\underline{\underline{\mathbf{B}}}\,\underline{\underline{\mathbf{B}}}^{\prime-1}\underline{\underline{\mathbf{A}}}^{\prime}\,\underline{\underline{\mathbf{A}}}^{-1}\right),$
$\displaystyle\underline{\underline{\mathbf{\varepsilon}}}_{p}$
$\displaystyle\approx$
$\displaystyle\frac{1}{2}\left(\left(\underline{\underline{\mathbf{B}}}\,\underline{\underline{\mathbf{B}}}^{\prime-1}\underline{\underline{\mathbf{A}}}^{\prime}\,\underline{\underline{\mathbf{A}}}^{-1}\right)^{T}\left(\underline{\underline{\mathbf{B}}}\,\underline{\underline{\mathbf{B}}}^{\prime-1}\underline{\underline{\mathbf{A}}}^{\prime}\,\underline{\underline{\mathbf{A}}}^{-1}\right)-\underline{\underline{\mathbf{1}}}\right),$

where the last line is a linear approximation. Alternatively we can express
Eq. LABEL:eqn:plasticStrainTensor in terms of the number of lattice repeats.

$\underline{\underline{\mathbf{\varepsilon}}}_{p}\approx\frac{1}{2}\left(\left(\underline{\underline{\mathbf{B}}}\,\underline{\underline{\mathbf{N}}}^{\prime}\,\underline{\underline{\mathbf{N}}}^{-1}\underline{\underline{\mathbf{B}}}^{-1}\right)+\left(\underline{\underline{\mathbf{B}}}\,\underline{\underline{\mathbf{N}}}^{\prime}\,\underline{\underline{\mathbf{N}}}^{-1}\underline{\underline{\mathbf{B}}}^{-1}\right)^{T}\right)-\underline{\underline{\mathbf{1}}}$
(6)

We conclude that the deformation corresponding to plastic strain is directly
related to the change in the reference lattice for the simulation, and can be
simply computed from the change in number of lattice repeats. This implies
that both transformations which shear the reference lattice, but preserve the
number of lattice sites
$(\mathrm{det}[\underline{\underline{\mathbf{N}}}^{\prime}]=\mathrm{det}[\underline{\underline{\mathbf{N}}}])$
_and_ transformations which change the number of lattice sites
$(\mathrm{det}[\underline{\underline{\mathbf{N}}}^{\prime}]\neq\mathrm{det}[\underline{\underline{\mathbf{N}}}])$
can be represented within the same plastic deformation
$\underline{\underline{\mathbf{F}}}_{p}$.

Hence in this work, the formalism for the plastic deformation gradient
$\underline{\underline{\mathbf{F}}}_{p}$ contains both radiation-induced slip
and swelling.

Now consider a general plastic transformation defined by a Burgers vector
change $\mathbf{b}$ each time we move along the line $\mathbf{n}$, with no
associated homogeneous lattice strain. The magnitude of the Burgers vector is
$b=|\mathbf{b}|$. We can define the degree of the transformation with a single
number $\gamma=b/|\mathbf{n}|$. $\gamma$ is simply the angle of shear if
$\mathbf{b}$ is in the plane with normal $\mathbf{n}$. Then the periodic
repeat vector $\mathbf{A}_{1}$ changes to
$\mathbf{A}_{1}+(\mathbf{A}_{1}\cdot\hat{\mathbf{n}}/|\mathbf{n}|)\mathbf{b}$,
where the caret $\hat{\mathbf{}}$, denotes a normal vector, and so the
supercell changes to

$\underline{\underline{\mathbf{A}}}^{\prime}=\underline{\underline{\mathbf{A}}}+\frac{\gamma}{b}\left(\begin{array}[]{ccc}(\mathbf{A}_{1}\cdot\hat{\mathbf{n}})\,\mathbf{b}&(\mathbf{A}_{2}\cdot\hat{\mathbf{n}})\,\mathbf{b}&(\mathbf{A}_{3}\cdot\hat{\mathbf{n}})\,\mathbf{b}\\\
|&|&|\\\ |&|&|\end{array}\right).$ (7)

It is convenient at this point to consider as an illustrative example a cubic
simulation cell, of side $ma_{0}$, i.e.
$\underline{\underline{\mathbf{A}}}=m\,a_{0}\underline{\underline{\mathbf{1}}}$.
Then the change in such a supercell due to a slip event is with Burgers vector
$\mathbf{b}$ on plane $\hat{\mathbf{n}}$ is given by the outer product

$\underline{\underline{\mathbf{A}}}^{\prime}=\underline{\underline{\mathbf{A}}}+\frac{m\gamma}{b}\mathbf{b}\otimes\hat{\mathbf{n}},$
(8)

and the linearised plastic strain is recognised as the Schmidt tensor,

$\underline{\underline{\mathbf{\varepsilon}}}_{p}=\frac{\gamma}{2b}\left(\mathbf{b}\otimes\hat{\mathbf{n}}+\hat{\mathbf{n}}\otimes\mathbf{b}\right).$
(9)

The new number of lattice repeats,
$\underline{\underline{\mathbf{N}}}^{\prime}=\underline{\underline{\mathbf{B}}}^{-1}\underline{\underline{\mathbf{A}}}^{\prime}=\underline{\underline{\mathbf{N}}}+(\gamma/b)\,\underline{\underline{\mathbf{N}}}\,(\mathbf{b}\otimes\hat{\mathbf{n}})$.

For BCC metals, $N_{\mathrm{motif}}=1$, so

$\displaystyle\underline{\underline{\mathbf{B}}}^{\mathrm{(BCC)}}$
$\displaystyle=$
$\displaystyle\frac{a_{0}}{2}\left(\begin{array}[]{ccc}-1&1&1\\\ 1&-1&1\\\
1&1&-1\\\ \end{array}\right),$ (13)
$\displaystyle\underline{\underline{\mathbf{N}}}^{\mathrm{(BCC)}}$
$\displaystyle=$ $\displaystyle\left(\begin{array}[]{ccc}0&m&m\\\ m&0&m\\\
m&m&0\\\ \end{array}\right),$ (17)

and the new number of lattice repeats in the BCC crystal due to a plastic
event is

${\underline{\underline{\mathbf{N}}}^{\mathrm{(BCC)}}}^{\prime}=\underline{\underline{\mathbf{N}}}^{\mathrm{(BCC)}}\\\
+\frac{m\gamma}{b}\left(\begin{array}[]{ccc}(b_{2}+b_{3})n_{1}&(b_{2}+b_{3})n_{2}&(b_{2}+b_{3})n_{3}\\\
(b_{3}+b_{1})n_{1}&(b_{3}+b_{1})n_{2}&(b_{3}+b_{1})n_{3}\\\
(b_{1}+b_{2})n_{1}&(b_{1}+b_{2})n_{2}&(b_{1}+b_{2})n_{3}\end{array}\right),$
(18)

while a similar process using an FCC primitive cell gives

${\underline{\underline{\mathbf{N}}}^{\mathrm{(FCC)}}}^{\prime}=\underline{\underline{\mathbf{N}}}^{\mathrm{(FCC)}}+\frac{m\gamma}{b}\left(\begin{array}[]{ccc}(-b_{1}+b_{2}+b_{3})n_{1}&(-b_{1}+b_{2}+b_{3})n_{2}&(-b_{1}+b_{2}+b_{3})n_{3}\\\
(b_{1}-b_{2}+b_{3})n_{1}&(b_{1}-b_{2}+b_{3})n_{2}&(b_{1}-b_{2}+b_{3})n_{3}\\\
(b_{1}+b_{2}-b_{3})n_{1}&(b_{1}+b_{2}-b_{3})n_{2}&(b_{1}+b_{2}-b_{3})n_{3}\end{array}\right)$
(19)

Note that, if we impose the condition that the reference lattice is a single
crystal, then the number of lattice repeats must be integer after a plastic
event, this constrains permissible values of $\gamma$ within the simulation
box. In particular, we see the smallest $\gamma$ must scale as $\sim 1/m$, so
small plastic deformations are not permitted in a small simulation box.

For the BCC/FCC cases, we find the change in lattice sites, $\Delta
N_{\mathrm{latt}}=\mathrm{Det}[\underline{\underline{\mathbf{N}}}^{\prime}]-\mathrm{Det}[\underline{\underline{\mathbf{N}}}]$,
is given by

$\displaystyle\Delta N_{\mathrm{latt}}^{\mathrm{(BCC)}}$ $\displaystyle=$
$\displaystyle 2\frac{m^{3}\gamma}{b}\mathbf{b}\cdot\hat{\mathbf{n}},$
$\displaystyle\Delta N_{\mathrm{latt}}^{\mathrm{(FCC)}}$ $\displaystyle=$
$\displaystyle 4\frac{m^{3}\gamma}{b}\mathbf{b}\cdot\hat{\mathbf{n}},$ (20)

from which we confirm that pure slip events
$(\mathbf{b}\cdot\hat{\mathbf{n}}=0)$ are associated with no change in lattice
sites, but other plastic events change the number of lattice sites and so
cause swelling. In Fig. 1, we illustrate the relationship between the slip
direction normal and the Burgers vector. Establishing a connection between
atomistic observations and continuum slip events offers the potential to
develop a comprehensive interpretative model for high-dose irradiation
scenarios. In section III.2.2, we outline both the similarities and
distinctions between atomic-level slip phenomena and macroscopic models of
plastic deformation.

Figure 1: Schematic representation depicting single crystal slip and
corresponding rigid rotation. (a) Prior to irradiation, (b) during
irradiation, as deformation initiates, (c) a pure slip event, and (d)
spontaneous atomic system rotation, within the constraint of a simulation box
disallowing shear.

###  II.3 Computing the elasto-plastic strain decomposition from atomic
positions

The reference primitive lattice vectors making the matrix
$\underline{\underline{\mathbf{B}}}$ can be found directly from the position
of peaks in the diffraction pattern. But this is an expensive calculation, so
below we describe an order-N approximation for
$\underline{\underline{\mathbf{B}}}$, which we find works well in the high-
dose irradiation cases described here.

Given a reference lattice $\underline{\underline{\mathbf{B}}}$, the reference
atom positions are

$\mathbf{y}=\underline{\underline{\mathbf{B}}}[u,v,w]^{T}+\mathbf{m}_{t},$
(21)

where $[u,v,w]$ is a triplet of integers and $\mathbf{m}_{t}$ is the position
of the $t^{th}$ motif point in the reference unit cell. Given the observed
atom positions, $\\{\mathbf{x}_{j}\\}$ , $j\in\\{1,\ldots N\\}$, the Wigner-
Seitz occupation, $o_{i}$, is the number of atoms closest to of reference atom
position $i$, i.e. the number of atoms for which
$\left|\mathbf{x}_{j}-\mathbf{y}_{i}\right|<\left|\mathbf{x}_{j}-\mathbf{y}_{k}\right|\,\forall\,k\neq
i$. The index $i\in\\{1,\ldots N_{\mathrm{latt}}\\}$, as defined above. Note
that the sum of the occupation, $\sum_{i}o_{i}=N_{\mathrm{atoms}}$, but that
$N_{\mathrm{atoms}}$ is not necessarily equal to $N_{\mathrm{latt}}$. The
number of point defects is defined by
$N_{\mathrm{PD}}=1/2\,\sum_{i}|o_{i}-1|$.

We can make the _ansatz_ that there exists a set of reference atom positions
$\\{\mathbf{y}_{i}\\}$ for which $N_{\mathrm{PD}}$ is minimised. For a single
crystal, we can in principle find this reference by minimising
$N_{\mathrm{PD}}$ with respect to a vector offset common to the motif points
$\\{\mathbf{m}_{t}\\}$ and the elements of
$\underline{\underline{\mathbf{N}}}$.

We can find a new reference lattice
$(\underline{\underline{\mathbf{B}}}^{\prime},\underline{\underline{\mathbf{N}}}^{\prime})$
from atom positions as follows. We compute the fit function on atom $i$:

$S_{i}=\underline{\underline{\mathbf{T}}}_{i}\sum_{j}\left(\mathbf{x}_{j}-\mathbf{x}_{i}-\delta\mathbf{y}_{k}\right)+\mathbf{\delta}_{i},$
(22)

where
$\delta\mathbf{y}=\underline{\underline{\mathbf{B}}}[u^{\prime},v^{\prime},w^{\prime}]^{T}+\mathbf{m}_{t}-\mathbf{m}_{1}$
is an expected separation between reference atom sites. This we minimise with
respect to the nine matrix elements of
$\underline{\underline{\mathbf{T}}}_{i}$, displacement vector
$\mathbf{\delta}_{i}$, and matching of neighbours to
$\\{\delta\mathbf{y}_{k}\\}$. We take the sum over neighbours with cutoff
between 3rd and 4th nearest neighbours. This method is slower than polyhedral
template matching [53], but robust when point defects are present. The matrix
$\underline{\underline{\mathbf{T}}}_{i}=\underline{\underline{\mathbf{R}}}_{i}(1+\underline{\underline{\mathbf{\varepsilon}}}_{i})$
is the deformation gradient, a product of atomic strain and rotation at atom
$i$. The rotation matrix and strain can be separated using polar
decomposition. To compute an approximate homogeneous lattice strain and
orientation, we need an appropriate average of the spatially varying
deformation gradient. To do this, we compute an arithmetic mean atomic
deformation gradient
$\langle\underline{\underline{\mathbf{T}}}\rangle=1/N\sum_{i}\underline{\underline{\mathbf{T}}}_{i}$,
and mean strain
$\langle\underline{\underline{\mathbf{\varepsilon}}}\rangle=1/N\sum_{i}\underline{\underline{\mathbf{\varepsilon}}}_{i}$
separately. Then we compute an appropriate rotation matrix,
$\underline{\underline{\mathbf{R}}}$, from the arithmetic mean using polar
decomposition. Finally we compute the homogeneous atomic deformation gradient
compatible with the unitary rotation matrix and mean strain, as
$\underline{\underline{\mathbf{T}}}=\underline{\underline{\mathbf{R}}}(1+\langle\underline{\underline{\mathbf{\varepsilon}}}\rangle)$.

The method above finds new reference lattice vectors
$\underline{\underline{\mathbf{B}}}^{\prime}=\underline{\underline{\mathbf{T}}}\,\underline{\underline{\mathbf{B}}}$
in $\mathcal{O}(N)$ time, and from this we find
$\underline{\underline{\mathbf{N}}}^{\prime}=\underline{\underline{\mathbf{A}}}^{\prime}\underline{\underline{\mathbf{B}}}^{\prime-1}$
to complete the reference lattice. We then search for an optimal vector offset
in $\\{\mathbf{m}_{t}\\}$. The Wigner-Seitz occupations $\\{o_{i}\\}$ reported
here therefore minimise $N_{\mathrm{PD}}$ subject to a $\mathcal{O}(N)$
estimation of the atomic deformation gradient.

###  II.4 Further analysis methods

The “Open Visualization Tool” (OVITO) [54] was used for visualization and the
dislocation extraction algorithm (DXA) [55] implemented was used for
dislocation identification. For the results obtained with our W-S analysis, a
cluster size distribution analysis was performed. In cases where a cluster
contains both vacancies and interstitials, we determined the net content of
defects to accurately assess the defect size. It is important to highlight
that such mixed clusters were infrequent in our study because of our process
of finding the reference lattice described above, and the disparity in the
cluster’s overall content before and after calculating the net defects was
approximately 1%, even at the higher doses. The cutoff value for vacancies and
interstitials are intermediate value between the nearest neighbor (1NN) and
2NN, and between the 3NN and 4NN, respectively [56, 57]. The shorter cutoff
for vacancies compared to previous studies was chosen to give an estimate for
possible deuterium retention, as for vacancies in clusters to collectively add
to the empty volume they need to be in the nearest position to another
vacancy.

##  III Results and Discussion

###  III.1 Defect evolution during irradiation

In this section, we highlight the core question of this paper by contrasting
the similarities and differences in the accumulation of radiation-induced
damage through various characterization methods. Fig. 2(a) shows the evolution
of point defects as a function of dose. The atomic fraction of vacancies,
computed from void isosurfaces using VID, is similar to the atomic fraction of
unoccupied lattice sites ($o_{i}=0$) computed with Wigner-Seitz below
$10^{-3}$ cdpa, but the levels diverge above this dose. Fig. 2(b) shows the
difference between these two methods. Unoccupied lattice sites which do not
correspond to an open volume are associated most commonly with vacancy loops.
This can be seen by comparing Fig. 2(c) and Fig. 2(d). The former, computed
with W-S, clearly shows vacancy loops bounded by dislocation lines in addition
to monovacancies and vacancy clusters, whereas the latter, computed with VID,
shows only monovacancies and vacancy clusters. Multiply occupied Wigner-Seitz
lattice sites ($o_{i}>1$) are associated with interstitials, either
indivdually as crowdions or in interstitial clusters and dislocation loops. We
see in Fig. 2(a) that the count of Wigner-Seitz unoccupied and multiply
occupied lattice sites diverges above $10^{-2}$ cdpa. Previous studies have
shown that similar apparent artefacts by W-S might appear if the perfect
lattice structure of a whole layer is identified as interstitials [58]. In
Fig. 2(c), the alignment of vacancies by W-S and dislocations by DXA confirms
that the significant increase in vacancy concentration observed in our work is
not an artefact. Hence, we pose the following question: What factors
contribute to the emergence of these vacancy loops and the divergence in the
evolution of vacancies and interstitials? This question will be answered in
detail in section III.2.2. Moreover, a full description of defect evolution
and clustering under our simulation conditions and a convergence study for the
CA can be found in the SI (see Table 1 and Figs. S1–S8).

Figure 2: (a) Point defect evolution as a function of dose characterized by
different methods. (b) The concentration of vacancies in loops calculated as
the difference between total unoccupied sites by the W-S method and vacancies
using the VID method. Solid lines: full MD simulations. Red squares: after
1600 CA steps. (c) Vacancy distribution obtained by W-S and (d) the void
distribution obtained by VID after CA, cdpa = 0.2047.

###  III.2 Microscopic deformation processes contributing to atomic and
plastic strain

In this section, we will illustrate swelling and dislocation creep phenomena
and their associated atomic and plastic strain responses that take place
during the irradiation process.

First we consider the average atomic strain, computed as
$1/3\,\mathrm{Tr}[\varepsilon_{a}]$. The evolution as a function of dose is
shown in Fig. 3. These results are qualitatively similar to the experimental
and simulation results reported in Ref. [50]. The atomic lattice strain is
tensile at low dose, during the accumulation of interstitial clusters and
vacancies. Interstitials have a large positive relaxation volume, while the
vacancies have a small negative relaxation volume, so if the numbers of each
balance, the strain is positive. At high dose, the interstitials have
coalesced into a network, leaving behind a vacancy-dominated microstructure
with compressive atomic strain. Here we show an improved quantitative match to
the experimental strain compared to Ref. [50]. This improvement is largely
because we use MD cascade annealing here, whereas the previous work used CRA
only. We should, however, note that the experimental boundary conditions were
different - the experimental study considered self-ion irradiation into a
slab, and so expansion was only possible in the $z$-direction. The micro-Laue
x-ray measurement detected $\varepsilon_{a}$ in the $z$-direction. As noted
above, in this work the simulation cell was allowed to expand in the $x$-,
$y$-, and $z$-directions. The change of sign in the strain at a dose of 0.01
cdpa corresponds to the point in Fig. 2(a) where the number of lattice sites
is changing, and the count of unoccupied and multiply-occupied W-S sites
diverges. In Ref. [50] this was attributed to interstitial plane formation and
here we consider the change in atomic and plastic deformation more generally.

The full evolution for both atomic and plastic deformation, rotation and the
other measurables are found in the SI, for all our investigated samples. In
the following subsections we are going through the main mechanisms and their
implication for selected examples.

Figure 3: Lattice strain derived from simulations. The gray shaded area
represents the error fluctuation range of the full MD section, showing the
range of three simulations. The error bars for the CRA+CA points is also
associated with the three different runs. The experimental strain data is from
Ref. [50].

####  III.2.1 Atomic deformation

We first consider the atomic strain in one of the simulations at a dose around
0.01 cdpa. At this point, the interstitial loops have grown large, but a
network has not formed. Fig. 4 shows the atomic strain $\varepsilon_{a}$ in
the $x$-, $y$-, and $z$-directions as a function of dose. In a small
simulation cell, it is quite expected that the components of the strain are
different, as they are determined by relatively few dislocation loops. In Fig.
4 we see the $x$\- and $z$-components are large below 0.01 cdpa, with the
$y$-component smaller. This is associated with the growth of a 1/2¡111¿
interstitial loop with habit plane (101), which can be seen in the top panel
of Fig. 4, i.e. snapshot A. Also visible in the snapshot is a smaller loop
with habit plane (-101). Over the course of 20 cascades (0.0001 cdpa), these
two loops join together to produce a mixed-habit-plane object seen in snapshot
B, and then coalesce to form a single large loop with habit plane (110) in
snapshot C. The dislocation line length is shown in the bottom panel of Fig.
4. There is no significant change in the fraction of W-S or VID defects during
this event (as can be seen from Figs. S2 and S3 in SI).

This habit plane rotation from (101) to (110) has a dramatic effect on the
_atomic_ strain components. In Fig. 4 middle panel we see the $y$-component
rapidly increase while the $z$-component decreases. However, this event is not
linked to a change in _plastic_ strain, as per the definition provided in this
paper. This lack of association stems from the unaltered matrix of lattice
repeats $\underline{\underline{\mathbf{N}}}$, as evidenced by the horizontal
dot line, coupled with the fact that all plastic components remain at a
constant value of 0. In the SI we provide several other examples of loop
coalescence and habit plane rotation leading to significant and rapid changes
in elastic stress, while not changing the plastic strain at all (see Figs.
S9–S15 and S17–S19).

These rapid fluctuations in local stress lead to rapid local changes in the
microstructure. It is unclear from the current simulations whether such rapid
local changes in the atomic stress will be so pronounced at greater length
scales comparable to grain sizes. These results therefore show the ease with
which the elastic dipole tensor of individual defects can be changed in
response to stress [59]. At low dose, the habit planes of defects can shift,
leading to a reorientation of their contribution to swelling [13].

Figure 4: Bottom row: the overall dislocation density evolution, with a
zoomed-in view of the period during which the strain drop occurs. Middle row
displays the evolution of atomic strain throughout the full MD simulation,
with a closer look at the strain drop during 20 rounds of irradiation. The
lattice sites have no change during the whole process, which indicates no
plastic transformation happens. To distinguish between three different states
of the evolution, a semi-transparent blue, yellow, and red background is used
in the zoomed-in figures. Upper row: to illustrate the strain transition
process, three representative frames are selected (A, B and C), namely before
the release of $\epsilon_{a}(z)$, during the release, and after the release.
The figures show the migration of interstitial loops at these three different
stages. Red arrows indicate the direction of loop movements, while red and
blue atoms represent vacancies and interstitials, respectively. Dislocation
lines with 1/2¡111¿ and ¡100¿ Burgers vectors are shown in green and pink
lines, respectively.

The key role of Eq. (20) lies in correlating atomistic observations with
continuum slip events, simplifying the comprehension that rotation of the
habit plane or a localized slip event does not impact the global cell,
maintaining a zero plastic response. Nonetheless, local slips occur in the
simulations, for instance, during the glide of dislocation lines/loops, as
depicted in Fig. 4. We formulated a straightforward slip tracing model to gain
deeper insights into the dislocation movement, emphasizing that this method is
an estimation, without meticulously converting each dislocation line network
into a continuous field that accurately signifies the network’s geometry and
Burgers vector alterations [60]. Notably, the minor slip amplitude from
cascade events allows for easy tracking of dislocation line movements without
employing the Sweep-tracing algorithm (STA), detailed in another study [60].
Fig. 5(a) schematically illustrates the slip tracing model. The core concept
of our model involves identifying the relevant slip facet by extracting
dislocation lines exhibiting glide. This is achieved by pairing segments one-
to-one before and after the slip event. The cross product of the slip vector
and the unit tangent vector of the dislocation line symbolizes the
$\hat{\mathbf{n}}_{i}$ for each slip facet. Fig. 5(b) displays the integral
statistics of slip facet normals calculated on each pair of dislocation loops
extracted from 20 full MD cascades, when the habit plane rotates from state A
to C in Fig. 4. The distribution of slip normals are plotted in projection as
a function of the azimutal ($\theta$) and elevation ($\phi$) angle, by color
mapping the density of the slip facet distributions. The high density slip
normals are identified as (11-2) and (01-1). In BCC metals, the planes with
the largest interplanar spacing are the 110 planes followed by the 112 and the
123 planes [61]. These slip planes often translates to lower energy barriers
for slip compared to other planes, especially under irradiation which can
enhance point defect mobility and dislocation line interaction. Other peaks
are also observed ((001)), which can be the contribution from dislocation
loops growing and rotating.

Figure 5: Slip behavior of dislocations during the full MD case. (a) Schematic
illustration of slip tracing extraction. The gray and black spheres symbolize
dislocation segments before and after a plastic slip event (with the interval
denoting a single cascade) extracted by the DXA method. The unit slip normal
is shown as $\hat{\mathbf{n}}$. (b) The cumulative distribution of
$\hat{\mathbf{n}}$ extracted over the course of 20 cascades from state A to C
in the full MD case.

####  III.2.2 Elasto-plastic deformation

To quantify the plastic strain, we proceed to examine the high-dose cases
through CRA+CA simulations, which result in highly deformed microstructures,
featuring abundant vacancies, offering the potential for plastic
transformation.

To begin with, the strain response at high doses is investigated. As an
example, we consider the 1st CA case with the highest starting dose of 0.2
cdpa. As shown in Fig. 6(a), the evolution of atomic strain $\varepsilon_{a}$
with both normal and shear components is observed. Initially, due to unstable
isolated defects introduced by the CRA, the cell exhibits an ill-defined
structure, leading to intense fluctuations in different indicators when the
dose is below 0.203 cdpa. Subsequently, the sudden increases or decreases in
excess lattice sites are seen, accompanied by abrupt changes in atomic strain
in specific directions. These corresponding changes in shear components
signify the occurrence of plastic transformation events, as shown in Fig.
6(b). The lattice strain changes instantly when the excess lattice sites form
or disappear, indicating that the plastic transformation happens immediately
after a single cascade event. We select two representative regions for a
detailed look, denoted as the blue and orange highlighted regions. These
particular sections are chosen due to abrupt observable alterations in atomic
strain. However, they exhibit different behavior: the region highlighted in
blue undergoes a more significant change, coinciding with a global plastic
slip event (Fig. 6(b)), which is not observed in the region highlighted in
orange.

Fig.6(c) and (d) show the atomic responses and cumulative slip normal
distributions for the two distinct cases, respectively. In Fig.6(c), the
instant change in atomic strain components is apparent, immediately following
a single cascade event, accompanying the onset of global plastic deformation.
An anomalous slip occurs in this instance (in addition to the conventional
123), as depicted by the cumulative distribution of slip normals in Fig.6(c).
Throughout this transformation, the non-conventional slip systems are
activated: the 111 slip normal prevails upon the occurrence of global plastic
deformation. The growth of dislocation networks provide additional slip
channels or influence the mobility and interaction of dislocations with 001
and 111 planes. Contrarily, Fig. 6(d) shows a 110-dominated slip duration when
there is no global plastic response but a change of atomic strain. Subsequent
examinations unveiled an additional instance of habit plane rotation. However,
the difference here compared to the low-dose case shown in Fig. 4 is that the
rotated dislocation loop is vacancy-type. Consequently, owing to the
relaxation volume of a vacancy loop is negative [62], vacancies distributed in
a material produce negative lattice strain, which can be readily observed
using x-ray diffraction [63, 64, 62, 50] in experiments. Therefore, an evident
opposite change occurs in the _atomic_ strain components compared to the
interstitial loop case shown in Fig. 4.

Analyzing the local slip event, it can be found that, even though the habit
plane rotation provides some additional slip channels, the slip planes should
still be dominated by the ones with largest interplanar spacings. When global
plastic deformation happens, this constraint no longer holds true. For the
region highlighted in blue, this process is associated with a change in the
number of lattice points, causing a plastic transformation characterized by
swelling events. Consequently, the count of interstitials and vacancies
diverges. This observed behavior is related to the Burgers vector $\mathbf{b}$
occurring across plane $n$, and it satisfies the relationship
$\mathbf{b}\cdot\hat{\mathbf{n}}\neq 0$, as indicated by Eq. (20). It is
important to highlight that this type of plastic event predominates in our
high-dose simulations, which might be a consequence of the orthorhombic
constraint imposed during the supercell relaxation. This observation suggests
that swelling tends to occur spontaneously when a plastic transformation is
triggered. The prevalence of this phenomenon underscores the strong
correlation between plasticity and the occurrence of swelling in our
simulations at high doses. The other cases experiencing plastic deformation
are shown in SI (see Figs. S16 and S20–S26).

Figure 6: (a) The elastic and (b) plastic responses as a function of dose at
the highest starting dose of 0.2 cdpa during the 1st CA. The amount of change
in lattice sites (according to the right axis) at different doses is indicated
by the red dashed line. The two representative change of atomic strain, i.e.
accompanying with global plastic deformation or not, are highlighted by the
semi-transparent blue and orange backgrounds, respectively. Closer looks at
local atomic strain and associated cumulative slip behavior during (c) global
plastic deformation (highlighted in blue) and (d) pure elastic response
(highlighted in orange).

Furthermore, our observations reveal the occurrence of rigid rotation during
plastic transformation. Notably, this rotation remains at zero during
instances of elastic deformation, such as in the full MD cases. However, here,
it is important to highlight that once the microstructure stabilizes,
typically at a dose value above 0.205 cdpa, rigid rotation persists even in
the absence of additional plastic slip events. Fig. 7 illustrates this
phenomenon, particularly in the zoomed-in detail as shown in Fig. 7(a), where
vacancy loop habit plane rotation induces symmetric evolution trends in the
$xy$\- and $xz$-shear components. Clearly, the $xz$ components exhibit a more
pronounced shear angle when compared to the $xy$ components. As a result, the
entire cell undergoes a slight rotation about the $y$-axis as the axis of
rotation, as depicted in Fig. 7(b) and (c). This phenomenon arises because
shear strain is not permitted within the LAMMPS simulation box in this study.
Consequently, when the orientation of the straining axis aligns favorably with
one of the crystal’s slip system, initiating slip along this primary slip
vector would distort the cell’s shape. However, the cell is constrained to
remain aligned with the simulation box, necessitating a rotational adjustment
to accommodate for this behavior. By monitoring the direction of rotation, we
can establish a connection with the coherence between continuum and atomic-
level slip events. Schmidt’s prediction suggests that, a crystal should
naturally rotate to align its predominant slip direction with the straining
axis, based on purely geometrical considerations [65, 66]. When comparing the
direction of rotation to that shown in Fig. 6(b), we observe a consistent
change in the direction of various shear components.

Figure 7: (a) Rigid rotation and amount of change in lattice sites (right axis
legend, red dashed line) at different doses for the highest starting dose (0.2
cdpa) during the 1st CA, $cos(\alpha_{ij})$ indicates the angle between the
shear component and the normal axes, where $i,j={x,y,z}$. (b) Corresponding
trajectories of rotation for [100] orientation, the doses are represented in a
cold to warm color scale. (c) Schematic illustration of the crystal rotation
(from the light grey to the black frame) of the atomic system with the dose
above 0.205 cdpa.

In the preceding discussion, we explored the plastic transformation in cases
of spontaneous swelling at high doses. Now, from Eq. (20), where we identify
another distinct microstructural event: when
$\mathbf{b}\cdot\hat{\mathbf{n}}=0$, indicating a pure slip event. Upon
careful examination, we discovered that such events are exceptionally rare.
For a pure slip phenomenon to occur, there must be a spontaneous change in
periodic lattice repeats, which only happens when the elements in the periodic
lattice repeat matrix undergo symmetrical position exchanges. During our
investigation, we encountered just one instance of this, at the very beginning
in the 3rd CA case, with the highest starting dose of 0.2 cdpa. As depicted in
Fig. 8(a), the light-blue-highlighted region in the plastic strain
$\varepsilon_{p}$ panel illustrates this pure slip event. Before and after
this single cascade event, both the normal and shear components exhibit
responses, but the number of excess lattice sites remains constant. A closer
examination reveals that the matrix of primitive cell repeats changes from

$\underline{\underline{\mathbf{N}}}_{\mathrm{before}}=\left(\begin{array}[]{ccc}0&121&121\\\
120&0&121\\\ 120&121&0\\\ \end{array}\right)\hskip 42.67912pt\mathrm{to}\hskip
42.67912pt\underline{\underline{\mathbf{N}}}_{\mathrm{after}}=\left(\begin{array}[]{ccc}0&121&121\\\
120&0&120\\\ 121&121&0\\\ \end{array}\right).$ (23)

To facilitate a direct comparison, we have annotated two cascade events,
namely labeled as regions A and B. The transition in region A corresponds to a
pure slip event, whereas the transition in region B represents a slip event
accompanied by lattice swelling. In order to gain a visual understanding of
the relationship between the angle defined by the Burgers vector and the slip
normals, dot product with the spatial Burgers vectors is performed to check
the perpendicularity. Fig. 8(b) shows the histogram of
$\mathbf{b}\cdot\hat{\mathbf{n}}$ after applying the Bootstrap method [67].
Employing normal distribution fitting, we obtained estimated result intervals
within the 95% confidence range for events shown in the region A and B
respectively. The central values for these two cases are approximately 0 and
7. Consequently, it is interesting to note that by discerning alterations in
the unit cell repeat number matrix, we can directly “observe” the occurrence
and type of slip events. Additionally, the infrequent pure slip events
observed in our simulations could be attributed to the boundary conditions
employed. It is conceivable that altering these conditions to include applied
shear stress during simulated irradiation might yield different results. In
our case, the practical result is that plastic events lead to material
swelling.

Figure 8: Direct observations of two distinct slip events. (a) Plastic strain
response (top panel) and corresponding defect concentration (bottom panel) as
a function of dose at the highest starting dose of 0.2 cdpa during the 3rd CA
case. The pure slip event and the slip event accompanied by lattice swelling
are highlighted by the semi-transparent blue and orange backgrounds,
respectively. The amount of change in lattice site (right axis legend) at
different doses is indicated by the red dashed line. (b) Schematic
illustration of dislocation slip. The gray and black spheres symbolize
dislocation segments before and after a plastic slip event (with the interval
denoting a single cascade) extracted by the DXA method. (c) Histogram of the
distribution of $\mathbf{b}\cdot\hat{\mathbf{n}}$ after the Bootstrap method
for events denoted as A and B. Normal distribution fitting curves are shown in
colors matching their respective data sets.

##  IV Conclusions

With the combination of full defect build-up and accelerated simulated
annealing MD simulations, we have investigated the mechanisms responsible for
atomic and plastic strain during damage accumulation and evolution. We devised
a strain decomposition method which works by resolving the homogeneous atomic
strain field from the atom positions, and using this to infer the plastic
strain required for a compatible single crystal. With our approach to
dissecting elasto-plastic deformation at the atomic level, we bring clarity to
the mechanism by which the gradual accumulation of lattice defects, as the
dose increases, leads to the development of a uniform lattice strain.

The main conclusions of this study are:

  1. 1.

At low dose, plastic strain changes are absent, but atomic strain changes are
brought about both by the accumulation of irradiation defects, and by the
response of the defects to the strain. This was illustrated by the mechanism
of habit plane rotation. In this work the defects were able to reduce elastic
strain energy in the simulation within molecular dynamics timescales,
indicating very low thermal activation barriers. In a real material,
impurities such as carbon may slow the movement of sections of the dislocation
loop [68]. It is unclear how the competition between the build-up of local
stress and the retarding effects of impurities are balanced in technologically
relevant conditions, and this merits further study.

  2. 2.

At high doses, both atomic and plastic strain can undergo alterations. The
plastic transformation can be categorized into two types: pure slip events and
slip accompanied by lattice swelling. The occurrence of pure slip events is
infrequent in our simulations. In cases where the number of lattice sites
remains unchanged, interstitials and vacancies show identical evolution
trends. However, when the plastic transformation related to swelling occurs, a
substantial change in excess lattice sites becomes apparent, which explains
why the interstitial and vacancy count diverge. In this case, anomalous slip
can occur. Lattice-swelling-related deformation could potentially activate
slip on these planes despite them not being the primary slip planes in BCC
tungsten, such as the 111 slip plane. The limited probability of pure slip
events in the studied dose range accounts for the consistent material swelling
observed during high-dose irradiation.

This study presents a straightforward and effective method for decoupling the
lattice strain field, which influences the evolution of defects resulting from
irradiation, from the plastic strain field, which leads to material swelling.
The methodology presented here can be used as an automatic analysis for future
simulations, which can identify the elasto-plastic deformation. The results of
this method can additionally be used to locate the interesting phenomena,
which enables analysis by other methods, which cannot be used globally or on
the the whole simulation series. This work is therefore a step towards
describing radiation defect accumulation within the framework of finite
element elasticity models. The findings underscore the detrimental impact of
high-dose irradiation, leading to significant creep phenomena that pose
challenges for the performance of advanced fission and fusion reactor
components.

## Acknowledgements

The authors would like to thank Max Boleininger and Luca Reali for helpful
discussions. This work has been carried out under the DEVHIS project, funded
by the Academy of Finland (Grant number 340538). This work has been partly
carried out within the framework of the EUROfusion Consortium, funded by the
European Union via the Euratom Research and Training Programme (Grant
Agreement No 101052200 — EUROfusion). Views and opinions expressed are however
those of the author(s) only and do not necessarily reflect those of the
European Union or the European Commission. Neither the European Union nor the
European Commission can be held responsible for them. Computer time granted by
the IT Center for Science – CSC – Finland and the Finnish Grid and Cloud
Infrastructure (persistent identifier urn:nbn:fi:research-infras-2016072533)
is gratefully acknowledged.

## References

  * ITER Physics Basis Editors [1999] ITER Physics Basis Editors, ITER Physics Basis, Nucl. Fusion 39, 2137 (1999).
  * Eckstein and Laszlo [1991] W. Eckstein and J. Laszlo, J. Nucl. Mater. 183, 19 (1991).
  * Lide [2005] D. R. Lide, ed.,  _CRC Handbook of Chemistry and Physics_ (CRC Press, Boca Raton, FL, 2005).
  * Wang  _et al._ [2011] L. G. Wang, A. Van De Walle, and D. Alfè, Phys. Rev. B 84, 092102 (2011).
  * Rieth  _et al._ [2011] M. Rieth, J. L. Boutard, S. L. Dudarev, T. Ahlgren, S. Antusch, N. Baluc, M. F. Barthe, C. S. Becquart, L. Ciupinski, J. B. Correia, C. Domain, J. Fikar, E. Fortuna, C. C. . Fu, E. Gaganidze, T. Galán, C. García-Rosales, B. Gludovatz, H. Greuner, K. Heinola, N. Holstein, N. Juslin, F. Koch, W. Krauss, K. Kurzydlowski, J. Linke, C. Linsmeier, N. Luzginova, H. Maier, M. Martínez, J. Missiaen, M. Muhammed, A. Muñoz, M. Muzyk, K. Nordlund, D. Nguyen-Manh, P. Norajitra, J. Opschoor, G. Pintsuk, R. Pippan, G. Ritz, L. Romaner, D. Rupp, R. Schäublin, J. Schlosser, I. Uytdenhouwen, J. van der Laan, L. Veleva, L. Ventelon, S. Wahlberg, F. Willaime, S. Wurster, and M. Yar, J. Nucl. Mater. 417, 463 (2011).
  * Garner  _et al._ [1981] F. A. Garner, M. L. Hamilton, N. F. Panayotou, and G. D. Johnson, J. Nucl. Mater. 104, 803 (1981).
  * Odette and Lucas [1998] G. R. Odette and G. E. Lucas, Radiat. Eff. Defects Solids 144, 189 (1998).
  * Zinkle  _et al._ [1993] S. J. Zinkle, P. J. Maziasz, and R. E. Stoller, J. Nucl. mater. 206, 266 (1993).
  * Monnet [2018] G. Monnet, J. Nucl. Mater. 508, 609 (2018).
  * Reza  _et al._ [2020] A. Reza, H. Yu, K. Mizohata, and F. Hofmann, Acta Mater. 193, 270 (2020).
  * Hollingsworth  _et al._ [2022] A. Hollingsworth, M. F. Barthe, M. Y. Lavrentiev, P. M. Derlet, S. L. Dudarev, D. R. Mason, Z. Hu, P. Desgardin, J. Hess, S. Davies, B. Thomas, H. Salter, E. F. J. Shelton, K. Heinola, K. Mizohata, A. De Backer, A. Baron-Wiechec, I. Jepu, Y. Zayachuk, A. Widdowson, E. Meslin, and A. Morellec, J. Nucl. Mater. 558, 153373 (2022).
  * Frydrych [2023] K. Frydrych, Crystals 13, 771 (2023).
  * Dudarev  _et al._ [2018] S. L. Dudarev, D. R. Mason, E. Tarleton, P. W. Ma, and A. E. Sand, Nucl. Fusion 58, 126002 (2018).
  * Xiao  _et al._ [2015] X. Xiao, D. Terentyev, L. Yu, D. Song, A. Bakaev, and H. Duan, J. Nucl, Mater. 466, 312 (2015).
  * Derlet and Dudarev [2020] P. M. Derlet and S. L. Dudarev, Phys. Rev. Mater. 4, 023605 (2020).
  * Boleininger  _et al._ [2023] M. Boleininger, D. R. Mason, A. E. Sand, and S. L. Dudarev, Sci Rep 13, 1684 (2023).
  * Mason  _et al._ [2021] D. R. Mason, F. Granberg, M. Boleininger, T. Schwarz-Selinger, K. Nordlund, and S. L. Dudarev, Phys. Rev. Mater. 5, 095403 (2021).
  * Granberg  _et al._ [2023] F. Granberg, D. R. Mason, and J. Byggmästar, Comput. Mater. Sci. 217, 111902 (2023).
  * Kiener  _et al._ [2011] D. Kiener, P. Hosemann, S. A. Maloy, and A. M. Minor, Nat. Mater. 10, 608 (2011).
  * Schiøtz  _et al._ [1999] J. Schiøtz, T. Vegge, F. D. Di Tolla, and K. W. Jacobsen, Phys. Rev. B 60, 11971 (1999).
  * Van Swygenhoven and Caro [1998] H. Van Swygenhoven and A. Caro, Phys. Rev. B 58, 11246 (1998).
  * Van Swygenhoven  _et al._ [1999] H. Van Swygenhoven, M. Spaczer, A. Caro, and D. Farkas, Phys. Rev. B 60, 22 (1999).
  * Van Swygenhoven  _et al._ [2001] H. Van Swygenhoven, A. Caro, and D. Farkas, Mater. Sci. Eng. A 309, 440 (2001).
  * Van Swygenhoven and Derlet [2001] H. Van Swygenhoven and P. M. Derlet, Phys. Rev. B 64, 224105 (2001).
  * Yamakov  _et al._ [2001] V. Yamakov, D. Wolf, M. Salazar, S. R. Phillpot, and H. Gleiter, Acta Mater. 49, 2713 (2001).
  * Yamakov  _et al._ [2002a] V. Yamakov, D. Wolf, S. R. Phillpot, and H. Gleiter, Acta Mater. 50, 5005 (2002a).
  * Yamakov  _et al._ [2002b] V. Yamakov, D. Wolf, S. R. Phillpot, A. K. Mukherjee, and H. Gleiter, Nat. Mater. 1, 45 (2002b).
  * Yamakov  _et al._ [2003] V. Yamakov, D. Wolf, S. R. Phillpot, and H. Gleiter, Acta Mater. 51, 4135 (2003).
  * Yamakov  _et al._ [2004] V. Yamakov, D. Wolf, S. R. Phillpot, A. K. Mukherjee, and H. Gleiter, Nat. Mater. 3, 43 (2004).
  * Vo  _et al._ [2008] N. Q. Vo, R. S. Averback, P. Bellon, S. Odunuga, and A. Caro, Phys. Rev. B 77, 134108 (2008).
  * Stukowski and Arsenlis [2012] A. Stukowski and A. Arsenlis, Modell. Simul. Mater. Sci. Eng. 20, 035012 (2012).
  * Nordlund  _et al._ [1998] K. Nordlund, M. Ghaly, R. S. Averback, M. Caturla, T. Diaz de La Rubia, and J. Tarus, Phys. Rev. B 57, 7556 (1998).
  * Leino  _et al._ [2018] A. A. Leino, G. D. Samolyuk, R. Sachan, F. Granberg, W. J. Weber, H. Bei, J. Liu, P. Zhai, and Y. Zhang, Acta Mater. 151, 191 (2018).
  * Bhardwaj  _et al._ [2021] U. Bhardwaj, A. E. Sand, and M. Warrier, Modell. Simul. Mater. Sci. Eng. 29, 065015 (2021).
  * Goryaeva  _et al._ [2020] A. M. Goryaeva, C. Lapointe, C. Dai, J. Dérès, J. B. Maillet, and M. C. Marinica, Nat. Commun. 11, 4691 (2020).
  * Anderson  _et al._ [2017] P. M. Anderson, J. P. Hirth, and J. Lothe,  _Theory of dislocations_ (Cambridge University Press, 2017).
  * Dudarev  _et al._ [2010] S. L. Dudarev, M. R. Gilbert, K. Arakawa, H. Mori, Z. Yao, M. L. Jenkins, and P. M. Derlet, Phys. Rev. B 81, 224107 (2010).
  * Mason  _et al._ [2014] D. R. Mason, X. Yi, M. A. Kirk, and S. L. Dudarev, J. Phys. Condens. Matter 26, 375701 (2014).
  * Plimpton [1995] S. Plimpton, J. Comput. Phys, 117, 1 (1995).
  * Ghaly  _et al._ [1999] M. Ghaly, K. Nordlund, and R. S. Averback, Phil. Mag. A 79, 795 (1999).
  * Nordlund [1995] K. Nordlund, Comp. Mater. Sci. 3, 448 (1995).
  * Ziegler [2004] J. F. Ziegler, Nucl. Instr. Meth. Phys. Res. Sec. B 219, 1027 (2004).
  * Sand  _et al._ [2013] A. E. Sand, S. L. Dudarev, and K. Nordlund, EPL (Europhysics Lett.) 103, 46003 (2013).
  * Ackland and Thetford [1987] G. J. Ackland and R. Thetford, Phil. Mag. A 56, 15 (1987).
  * Zhong  _et al._ [1998] Y. Zhong, K. Nordlund, M. Ghaly, and R. S. Averback, Phys. Rev. B 58, 2361 (1998).
  * Granberg  _et al._ [2021] F. Granberg, J. Byggmästar, and K. Nordlund, J. Nucl. Mater. 556, 153158 (2021).
  * Kinchin and Pease [1955] G. H. Kinchin and R. S. Pease, Reports on progress in physics 18, 1 (1955).
  * Robinson and Torrens [1974] M. T. Robinson and I. M. Torrens, Phys. Rev. B 9, 5008 (1974).
  * Norgett  _et al._ [1975] M. J. Norgett, M. T. Robinson, and I. M. Torrens, Nucl. Eng. Des. 33, 50 (1975).
  * Mason  _et al._ [2020] D. R. Mason, S. Das, P. M. Derlet, S. L. Dudarev, A. J. London, H. Yu, N. W. Phillips, D. Yang, K. Mizohata, R. Xu, and F. Hofman, Phys. Rev. Lett. 125, 225503 (2020).
  * Mura [2013] T. Mura,  _Micromechanics of defects in solids_ (Springer Science & Business Media, 2013).
  * Note [1] To see that the primitive cell can be rotated, we use proof by _reductio ad absurdum_ : The most extreme irradiation disordering could be taking out all the atoms, replacing them in amorphous positions, and annealing back to a crystal lattice. For no rotation to be true it is necessary that the exact same lattice orientation is always recovered from amorphisation + annealing.
  * Larsen  _et al._ [2016] P. M. Larsen, S. Schmidt, and J. Schiøtz, Modell. Simul. Mater. Sci. Eng. 24, 055007 (2016).
  * Stukowski [2009] A. Stukowski, Modell. Simul. Mater. Sci. Eng. 18, 015012 (2009).
  * Stukowski  _et al._ [2012] A. Stukowski, V. V. Bulatov, and A. Arsenlis, Modell. Simul. Mater. Sci. Eng. 20, 085007 (2012).
  * Terentyev  _et al._ [2006] D. Terentyev, C. Lagerstedt, P. Olsson, K. Nordlund, J. Wallenius, C. S. Becquart, and L. Malerba, J. Nucl. Mater. 351, 65 (2006).
  * Björkas and Nordlund [2007] C. Björkas and K. Nordlund, Nucl. Instr. Meth. Phys. Res. Sec. B 259, 853 (2007).
  * Wu  _et al._ [2021] J. Wu, Z. Xu, L. Liu, A. Hartmaier, M. Rommel, K. Nordlund, T. Wang, R. Janisch, and J. Zhao, J. Mat. Chem. C 9, 2258 (2021).
  * Dudarev and Sutton [2017] S. L. Dudarev and A. P. Sutton, Acta Mater. 125, 425 (2017).
  * Bertin  _et al._ [2022] N. Bertin, L. A. Zepeda-Ruiz, and V. V. Bulatov, Mater. Theory 6, 1 (2022).
  * Weinberger  _et al._ [2013] C. R. Weinberger, B. L. Boyce, and C. C. Battaile, Int. Mater. Rev. 58, 296 (2013).
  * Mason  _et al._ [2019] D. R. Mason, D. Nguyen-Manh, M. C. Marinica, R. Alexander, A. E. Sand, and . L. Dudarev, J. Appl. Phys. 126 (2019).
  * Simmons and Balluffi [1960] R. O. Simmons and R. W. Balluffi, Phys. Rev. 117, 52 (1960).
  * Hertz  _et al._ [1973] W. Hertz, W. Waidelich, and H. Peisl, Phys. Lett. A 43, 289 (1973).
  * Meyers and Chawla [2008] M. A. Meyers and K. K. Chawla,  _Mechanical behavior of materials_ (Cambridge university press, 2008).
  * Zepeda-Ruiz  _et al._ [2021] L. A. Zepeda-Ruiz, A. Stukowski, T. Oppelstrup, N. Bertin, N. R. Barton, R. Freitas, and V. V. Bulatov, Nat. Mater. 20, 315 (2021).
  * Efron [1992] B. Efron, in  _Breakthroughs in statistics: Methodology and distribution_ (Springer, 1992) pp. 569–593.
  * Castin  _et al._ [2019] N. Castin, A. Dubinko, G. Bonny, A. Bakaev, J. Likonen, A. De Backer, A. E. Sand, K. Heinola, and D. Terentyev, J. Nucl. Mater. 527, 151808 (2019).