[1]

inst1]organization=UK Atomic Energy Authority,addressline=Culham Science
Centre, city=Abingdon, postcode=OX14 3DB, country=UK

inst2]organization=Department of Materials, University of
Manchester,city=Manchester, postcode=M13 9PL, country=UK

inst5]organization=Department of Materials Physics, Eötvös
University,postcode=PO Box 32, H-1518, city=Budapest, country=Hungary

inst6]organization=Deutsches Elektronen-Synchrotron DESY,addressline =
Notkestr. 85, postcode=22607, city=Hamburg, country=Germany

inst7]organization=Monash University, city=Clayton, postcode =VIC 3800,
country = Australia

[1]Corresponding author

# Dislocation density transients and saturation in irradiated zirconium

Andrew R. Warwick<EMAIL_ADDRESS>Rhys Thomas  M. Boleininger  Ö. Koç
G. Zilahi  G. Ribárik  Z. Hegedues  U. Lienert  T. Ungar  C. Race  M. Preuss
P. Frankel  S. L. Dudarev  [  [  [  [  [

###### Abstract

Zirconium alloys are widely used as the fuel cladding material in pressurised
water reactors, accumulating a significant population of defects and
dislocations from exposure to neutrons. We present and interpret synchrotron
microbeam X-ray diffraction measurements of proton-irradiated Zircaloy-4,
where we identify a transient peak and the subsequent saturation of
dislocation density as a function of exposure. This is explained by direct
atomistic simulations showing that the observed variation of dislocation
density as a function of dose is a natural result of the evolution of the
dense defect and dislocation microstructure driven by the concurrent
generation of defects and their subsequent stress-driven relaxation. In the
dynamic equilibrium state of the material developing in the high dose limit,
the defect content distribution of the population of dislocation loops,
coexisting with the dislocation network, follows a power law with exponent
$\alpha\approx 2.2$. This corresponds to the power law exponent of
$\beta\approx 3.4$ for the distribution of loops as a function of their
diameter that compares favourably with the experimentally measured values of
$\beta$ in the range $3\leq\beta\leq 4$.

###### keywords:

zirconium irradiation dislocations defects

Dislocation density evolution in irradiated Zr has been predicted and
measured.

A transient peak and subsequent saturation in dislocation density has been
observed.

Dislocation loop diameters in heavily irradiated Zr are power law distributed.

##  1 Introduction

In the core of modern boiling (BWR) or pressurized (PWR) water reactors, the
uranium dioxide fuel assemblies are immersed in circulating pressurised water
and thus it is critical that only the heat produced by the fission reactions
is transported by the coolant and there is no contamination of the coolant
from the radioactive fuel itself. Hence, the fuel is cladded to protect the
reactor environment from contamination, be that during reactor operation or in
transit. From the design choices that date back over fifty years (Rickover et
al., 1975), zirconium alloys are currently employed as the uranium dioxide
fuel cladding in water-cooled reactors. Containing more than 95 wt% Zr, these
alloys are mostly pure zirconium, chosen for its low neutron absorption cross
section (Pomerance, 1951). Small amounts of Sn, Nb, Fe and/or Cr in the alloys
help protect against corrosion and improve structural integrity (Lemaignan,
2012; Onimus and Bechade, 2012).

Figure 1: Dislocation density as a function of dose observed experimentally
using microbeam synchrotron X-ray diffraction (XRD) measurements of proton-
irradiated Zircaloy-4, compared to predictions derived from simulations of
pure $\alpha$-Zr performed using the creation relaxation algorithm technique.
Four experimental samples with nominal doses of 0.1, 0.5, 1.0 and 2.0 dpa were
scanned across the variable dose regime. The simulation results are scaled by
a factor of 0.1 and averaged over three different interatomic potentials (see
text).

The elastic neutron scattering cross-section for a Zr nucleus is, however,
similar to that of other elements in the Periodic table (Sears, 2006), and a
prolonged exposure to neutron irradiation results in the accumulation of a
considerable amount of microscopic radiation defects generated from the atomic
recoils initiated by collisions with neutrons. This gives rise to the
deterioration of mechanical and physical properties and stimulates dimensional
changes (Holt, 1988; Onimus and Bechade, 2012). The high energy
$>1\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ neutrons produced by the fissile
uranium oxide fuel (Nicodemus and Staub, 1953) collide with Zr nuclei,
initiating collision cascades that displace atoms from their lattice sites,
rearrange the crystal structure, and generate crystal lattice defects (Domain
and Legris, 2005). Defects accumulate with increasing exposure to neutron
irradiation, in the form of pairs of self-interstitial and vacancy defects
(Frenkel pairs) as well as in the form of clusters of defects that eventually
coalesce into large-scale defects, such as dislocation loops and dislocations
(Warwick et al., 2021). At temperatures above
300-$350\text{\,}\mathrm{\SIUnitSymbolCelsius}$ where point defect diffusion
occurs at an appreciable rate comparable or faster rates than the dose rate,
it is important to consider random Brownian motion of defects to interfaces,
grain boundaries, dislocations or other point defect clusters, whereas at
lower temperatures the microstructure of accumulating defects is dominated by
other factors. Given the significance of the life-limiting effect of
structural changes on the properties of zirconium cladding, there is
significant research effort aimed at improving the performance of structural
zirconium alloys in the operating environment of a fission reactor (Adamson et
al., 2019; Zinkle and Was, 2013).

The exposure of a material to energetic particles is often quantified by the
notion of dose, typically expressed in the units of ‘displacement per atom’
(dpa). Dpa is a simple measure of exposure of a material to radiation and
represents the average number of times each atom in the material has been a
part of a Frenkel self-interstitial-vacancy defect pair. Typically, a
zirconium alloy cladding is exposed to $\sim
15\text{\,}\mathrm{d}\mathrm{p}\mathrm{a}$ over the five years of service
(Zinkle and Was, 2013). Producing an accurate safety case involves the
identification of the microstructure formed at a given dose, temperature, and
externally applied load.

In particular, the formation of dislocation loops and dislocations is known to
play a critical role in the resulting deterioration of the cladding’s
structural and mechanical properties. For instance, in the large dose limit,
high densities of defects and dislocations accumulate in the cladding, causing
embrittlement. Another important degradation mode is the so-called
‘irradiation-induced growth’ (IIG) that arises from the anisotropy of the
hexagonal close packed (hcp) crystal structure of zirconium ($\alpha$-Zr) and
zirconium alloys (Griffiths, 2020; Onimus et al., 2022). This crystal
structure is stable up to and beyond the reactor core’s operating temperature
range of $280\text{\,}\mathrm{\SIUnitSymbolCelsius}$ to
$350\text{\,}\mathrm{\SIUnitSymbolCelsius}$, exhibiting an hcp-bcc instability
only at significantly higher temperatures above
$860\text{\,}\mathrm{\SIUnitSymbolCelsius}$ (Willaime and Massobrio, 1989).

With increasing dose, the interstitial and vacancy-type dislocation loops with
the $\frac{1}{3}\langle 2\bar{1}\bar{1}0\rangle$ Burgers vectors, inhabiting
the prismatic crystallographic planes, form a population of the so-called ‘a
loops’. At relatively low temperatures, this strongly correlates with observed
elongation along the ‘a’ $\langle 2\bar{1}\bar{1}0\rangle$ and contraction
along the ‘c’ $\langle 0001\rangle$ crystallographic directions, saturating at
doses larger than $1\text{\,}\mathrm{d}\mathrm{p}\mathrm{a}$ (Holt, 1988). At
temperatures above $\sim 300\text{\,}\mathrm{\SIUnitSymbolCelsius}$ and doses
exceeding $\sim 5\text{\,}\mathrm{d}\mathrm{p}\mathrm{a}$, large vacancy-type
‘c loops’ with Burgers vectors $\frac{1}{2}\langle 0001\rangle$ or
$\frac{1}{6}\langle 20\bar{2}3\rangle$ appear in the basal crystallographic
planes. Accompanying this onset of formation of the c-type vacancy loops, the
magnitudes of a and c strains increase linearly with dose in a phenomenon
called ‘the breakaway growth’ (Choi and Kim, 2013). Whilst it is known that
dislocations have substantial elastic relaxation volume (Dudarev and Ma, 2018;
Boleininger et al., 2022) giving rise to significant dimensional changes
(Christensen et al., 2020), the present understanding of this IIG is mostly
phenomenological and existing models are unable to predict, from first
principles, the variation of the dislocation content consistent with, or at
least verified by, the experimental data.

Below we present experimental observations and predictions derived from
simulations, showing that the density of dislocations saturates in the
temperature and dose range where dimensional changes exhibit saturation. This
is summarised in Figure 1 illustrating the dislocation densities
experimentally measured in proton-irradiated Zircaloy-4
(Zr-1.5%Sn-0.2%Fe-0.1%Cr) together with the simulation data for irradiated
$\alpha$-Zr plotted as a function of dose. In agreement with observations
performed using ion irradiation (Yu et al., 2017), we find that the density of
dislocations evolves through a transient at doses
$<1\text{\,}\mathrm{d}\mathrm{p}\mathrm{a}$ before saturating at larger doses.
The microstructures produced by our simulations indicate that at very low
doses $\ll 1\text{\,}\mathrm{d}\mathrm{p}\mathrm{a}$, small dislocation loops
form, which subsequently grow and coalesce into a complex interconnected
dislocation network developing at moderate doses
$<1\text{\,}\mathrm{d}\mathrm{p}\mathrm{a}$. The dislocation network
eventually forms complete crystallographic planes (Mason et al., 2020;
Boleininger et al., 2023) and partly dissociates into large dislocation loops
at high doses $\gg 1\text{\,}\mathrm{d}\mathrm{p}\mathrm{a}$. The formation of
dislocation loops and dislocations is principally driven by the stress-
mediated clustering of self-interstitial atoms, where in the high dose limit
the self-interstitial cluster population size distribution follows a power law
$p(N)\propto 1/N^{2.2}$ developing in the saturation regime corresponding to
high neutron exposure. Below, we show that this fact, also confirmed by
experimental observations (Ungár et al., 2021), has important implications for
the interpretation of experimental observations of dislocation loops and for
our understanding of the dynamics of microstructural evolution in the
irradiated cladding.

Our manuscript is structured by detailing our methods in §2 and presenting our
results in §3 before summarising key conclusions in §4. The concept and
relevant definitions of dislocation density are discussed in §2.1. Details of
our experimental set-up and and an overview of the cmwp line profile analysis
software are given in §2.2 followed by an explanation of our simulation
method, its range of validity and our choice of settings in §2.3. We show how
the simulated microstructure correlates to the dislocation density profile
shown in Figure 1 in §3.1 in addition to characterising the evolution of the
power law distribution of dislocation sizes. Finally, the stored energy
associated with dislocation content that drives the changing microstructure as
a function of dose is investigated in §3.2.

##  2 Theory and Methodology

###  2.1 Dislocation density

There are multiple ways to define the density of dislocation lines $\rho$ in a
deformed or irradiated material. For all the measurements and computations in
this study, $\rho$ is defined as a scalar ratio of the total dislocation line
length to volume $V$ containing the dislocations, namely (Hull and Bacon,
2011)

$\displaystyle\rho:=\frac{1}{V}\int_{\perp\in
V}{\left|\mathrm{d}\mathbf{l}\right|}.$ (1)

Here $\mathbf{l}$ is a position vector on a dislocation line such that
$\mathrm{d}\mathbf{l}=\bm{\xi}\mathrm{d}s$ for unit tangent vector $\bm{\xi}$
and arc length $s$, and the integration is performed with respect to the arc
length over all the dislocation lines $\perp$ in $V$. The choice of volume is
somewhat arbitrary but, for a given resolution, is expected to reflect the
average amount of dislocations present. For example, in a TEM micrograph this
can be chosen to be a region contained in the image and in a molecular
dynamics simulation one may use the entire simulation cell.

Another possible definition is the areal density $\rho_{A}$ that measures the
number of dislocation lines crossing an open surface as (Hull and Bacon, 2011)

$\displaystyle\rho_{A}:=\frac{1}{A}\int_{\perp\in
A}{\left|\mathrm{d}\mathbf{S}\cdot\bm{\xi}\right|},$ (2)

where $\mathrm{d}\mathbf{S}$ is a vector area element of $A$ with direction
normal to the surface and, similar to Equation 1, there is a dependence on the
choice of the surface. If all the dislocations are co-linear and perpendicular
to the chosen surface, Equations 1 and 2 produce the same value. Whilst for an
arbitrary distribution of dislocations this is often not the case, generally
both measures do not differ by much more than a factor of two (Schoeck, 1962).

Needless to say, a dislocation is not solely defined by its tangent vector and
it is noteworthy that neither Equation 1 nor Equation 2 contain any
information about the Burgers vector $\mathbf{b}$ of the dislocation.
Significant physical quantities such as dislocation energy and the Peach-
Koehler force both depend on $\mathbf{b}$. The Nye tensor $\bm{\alpha}$ is a
tensorial measure of dislocation density that is a linear function of lattice
curvature (Nye, 1953) and is also a function of position $\mathbf{x}$ such
that (Jones et al., 2016)

$\displaystyle\bm{\alpha}(\mathbf{x}):=\int_{\perp}{\delta\left(\mathbf{x}-\mathbf{l}\right)\mathbf{b}(\mathbf{l})\otimes\mathrm{d}\mathbf{l}},$
(3)

where $\otimes$ denotes the tensor product, integration is performed over all
of the dislocation lines in the system, and $\delta\left(\mathbf{x}\right)$ is
the Dirac delta distribution defined by the property

$\displaystyle\int\mathrm{d}^{3}\mathbf{x}^{\prime}\,\delta(\mathbf{x}^{\prime}-\mathbf{x})f(\mathbf{x}^{\prime})=f(\mathbf{x}),$
(4)

for an arbitrary well-behaved function $f(\mathbf{x})$. Whilst full
information about the dislocation content is contained in Equation 3,
attempting to average $\bm{\alpha}$ over a volume $V$ can be problematic.
Essentially, this stems from the fact that the integral of
$\mathrm{d}\mathbf{l}$ along a dislocation segment contained in $V$ that
starts and ends at $\mathbf{x}_{a}$ and $\mathbf{x}_{b}$ respectively is
$\mathbf{x}_{b}-\mathbf{x}_{a}$. Thus, any information pertaining to curvature
of a dislocation line is lost and closed paths in particular, i.e. dislocation
loops, provide no contribution to the volume average (Arsenlis and Parks,
1999; Mandadapu et al., 2014). The dislocations sections that integrate to
zero are referred to as Statistically Stored Dislocations (SSD) and the
surviving contributions are the Geometrically Necessary Dislocations (GND).
Experimental techniques that infer GND content, such as micro-beam Laue
measurements and high resolution electron back-scattered diffraction,
implicitly make use of Equation 3 (Das et al., 2018). As it is difficult to
characterise the entire population of dislocations using Equation 3,
throughout this manuscript we have chosen to use the scalar measure given by
$\rho$ in Equation 1 which is the same definition as that employed in our
X-ray line profile analysis.

###  2.2 Experiment

Dislocation densities were measured in four $3\text{\,}\mathrm{mm}$ $\times$
$1\text{\,}\mathrm{mm}$ $\times$ $0.5\text{\,}\mathrm{mm}$ Zircaloy-4 samples
(composition Zr-0.17Fe-1.24Sn-0.10Cr) proton-irradiated to different doses.
The samples possessed a recrystallised equiaxed microstructure with a low
dislocation density of $<1\text{\times}{10}^{14}\text{\,}\mathrm{m}^{-2}$ and
a characteristic ‘split-basal’ texture due to processing, where the basal
poles are aligned along the normal direction (ND), with a $\pm$30 degree tilt
towards the transverse direction. Irradiation induced growth strains in
similarly textured Zircaloys are known to saturate at doses below $\sim
10\text{\,}\mathrm{d}\mathrm{p}\mathrm{a}$ at
$320\text{\,}\mathrm{\SIUnitSymbolCelsius}$ (Adamson et al., 2019) and thus we
may expect dislocation densities to display a similar pattern of evolution in
these samples.

The ND face of each sample was proton irradiated with
$2\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ protons at
$350\text{\,}\mathrm{\SIUnitSymbolCelsius}$ at the University of Manchester’s
Dalton Cumbrian Facility, UK. The temperature of the samples during
irradiation was monitored via a thermal imaging camera in order to hold it
within $\pm 10\text{\,}\mathrm{\SIUnitSymbolCelsius}$ of the target
temperature. Unlike neutrons, the Coulomb interaction between the protons and
the target material results in the shallow penetration of protons into the
material. The resulting radiation exposure, quantified by the dose and dose
rate, varies significantly as a function of depth, with the dose rate being of
the order of $10^{-5}$ dpa/s. The dose profile was calculated using the quick
Kinchin-Pease setting in srim (Ziegler et al., 2010) with the lattice binding
energy and threshold displacement energy set to
$0\text{\,}\mathrm{e}\mathrm{V}$ and $40\text{\,}\mathrm{e}\mathrm{V}$
respectively (Stoller et al., 2013). A typical dose vs. depth profile in one
of our Zircaloy-4 samples consists of a plateau region extending $\sim
10\text{\,}\mathrm{\SIUnitSymbolMicro m}$ from the surface where the dose and
dose rate are approximately constant before sharply rising and falling to zero
at a region corresponding to protons coming to rest in the material, called
the Bragg peak and located at $\sim
30\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$. The samples were irradiated
such that the doses at 60% of the Bragg peak depth from the surface, termed
‘nominal doses’, were 0.1, 0.5, 1 and
$2\text{\,}\mathrm{d}\mathrm{p}\mathrm{a}$ respectively. Within the first
$30\text{\,}\mathrm{\SIUnitSymbolMicro m}$ of each sample from the irradiated
surface, the calculated dose and dose rate vary from their nominal values by
factors of 0.6 to 8.5 thus allowing us to measure data spanning over a wide
range of irradiation exposures.

Using a small X-ray beam of $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ $\times$
$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$ cross section, the samples were
scanned in cross-section from the surface to a depth of
$50\text{\,}\mathrm{\SIUnitSymbolMicro m}$ within the sample at
$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ increments at the P21.2 beamline at
the PETRA III synchrotron facility at DESY in Hamburg, Germany. The samples
were translated perpendicular to the scanning direction by
$200\text{\,}\mathrm{\SIUnitSymbolMicro m}$ during each scan to improve grain
statistics and reduce the spottiness of the pattern. The set-up of the
diffraction experiment and sample geometry is shown in Figure 2.

(a)

(b)

Figure 2:  (a) : Experimental geometry for measuring line profiles from
Zircaloy-4 samples on the P21.2 beamline at the DESY synchrotron, Hamburg,
Germany. (b): Variation of dislocation density and dose as a function of depth
for the nominal $0.1\text{\,}\mathrm{d}\mathrm{p}\mathrm{a}$ sample.

Dislocation densities were extracted from the line profiles using the
Convolutional Multiple Whole Profile (cmwp) software (Ribárik et al., 2020).
The cmwp software models the line profile intensity $I(q)$, where $q$ is the
wavevector magnitude, as a convolution of intensities arising from
instrumental effects, size broadening, and, in particular, strain broadening
due to dislocations. Instrumental effects were determined using a LaB6
standard specimen, and the size broadening was determined assuming a log-
normal size distribution of coherently scattering grains (Ribárik et al.,
2020). Here, we outline the model underpinning the cmwp software, whereas for
a broader context we refer an interested reader to reviews of the method
(Wilkens, 1970; Ungár et al., 1999; Ribárik et al., 2020).

In the theory of X-ray line profile analysis, the Fourier components of the
broadened intensity peak profile corresponding to a reciprocal lattice vector
$G$, denoted $I^{D}_{G}$, are related to the strain distribution by

$\displaystyle\mathcal{F}\left\\{I^{D}_{G}(q)\right\\}(L)=\exp\left[-2\pi^{2}G^{2}L^{2}\langle\epsilon^{2}_{G,L}\rangle\right],$
(5)

where $L$ is the Fourier variable and $\langle\epsilon^{2}_{G,L}\rangle$ is
the mean-square strain. Wilkens (1970) derived an expression for
$\langle\epsilon^{2}_{G,L}\rangle$ by numerical methods arising from the
‘restrictedly random distribution’ concept of dislocations. Dislocation lines
are assumed to be parallel in the sub-areas of equal size $A$ perpendicular to
the line direction. A number $N_{\perp}$ of dislocation lines with equal
numbers of positive and negative Burgers vectors occupy random positions in a
plane normal to the dislocation lines. The characteristic linear size of the
sub-areas is chosen to be proportional to a parameter termed the effective
outer cut-off radius $R_{e}$. The dislocation density $\rho$ is then defined
by Wilkens (1970) as the areal density of dislocations given by Equation 2
that, as discussed in §2.1, is equivalent to the volume density of
dislocations defined by Equation 1 for this specific dislocation
configuration. The dipole character of the distribution is determined by the
arrangement parameter $M=R_{e}\sqrt{\rho}$ . Whilst $R_{e}$ is one of the
fitting parameters in the cmwp software, thus affecting the value of $M$, this
article is principally concerned with the determination of $\rho$ and thus we
do not discuss the arrangement parameter further.

An expression for the mean-square strain for a restricted random distribution
of dislocations was derived by Wilkens (1970)

$\displaystyle\langle\epsilon^{2}_{G,L}\rangle=\frac{\rho
Cb^{2}}{4\pi}f(\eta),$ (6)

where $C$ is a parameter that is refined by the profile fitting cmwp
algorithm, termed the dislocation contrast factor. The value of $C$ was
evaluated for dislocations in $\alpha$-Zr by Balogh et al. (2016). The Wilkens
function $f(\eta)$, where $\eta=L/R_{e}$, has the following asymptotic forms
in the limit of small and large $\eta$, respectively

$\displaystyle f(\eta)\sim\begin{cases}\ln{\eta}&,\,\eta\rightarrow 0\\\
\frac{1}{\eta}&,\,\eta\rightarrow\infty.\end{cases}$ (7)

The numerically obtained formula for $f(\eta)$ may be found in Wilkens (1970).

Although the Wilkens model is formally derived assuming that a dislocation
configuration is composed of straight lines, the model is in fact able to
mimic the statistical properties of distributions of curved dislocations
(Kamminga and Delhez, 2000; Groma and Borbély, 2004). When compared to
transmission electron microscope (TEM) measurements of irradiated Zircaloy-2,
the cmwp software was able accurately follow the dislocation density evolution
as a function of dose (Seymour et al., 2017), and the cmwp approach has now
become an accepted tool for determining dislocation densities (Ungár et al.,
2021; Topping et al., 2018) as well as other microstructural features (Ungár
et al., 2021) in irradiated Zircaloys. The cmwp software evaluates parameters
describing effects of both the specimen size and dislocation broadening of
diffraction intensity peaks by first employing a statistical Monte Carlo
optimisation followed by the Marquardt-Levenberg non-linear least squares
algorithm, see Ribárik et al. (2020). All the peaks in the interval from
$0\text{\,}\mathrm{n}\mathrm{m}^{-1}$ $<q<$
$13\text{\,}\mathrm{n}\mathrm{m}^{-1}$ were included in the fitting procedure
and the uncertainty in $\rho$ was quantified according to the quality of fit
as described by Ribárik et al. (2020). The variation of dislocation density
with depth calculated by cmwp is shown in Figure 2(b) that, when mapped to the
dose calculated by srim at a given depth, enables plotting the dislocation
density as a function of dose as illustrated in Figure 1.

###  2.3 Simulation

The accumulation of defects and microstructural evolution as a function of
dose was simulated using the Creation Relaxation Algorithm (CRA) (Derlet and
Dudarev, 2020). The CRA exploits the separation of timescales associated with
relatively fast stress driven and comparatively slow thermally activated
evolution of defect microstructure. This results in a simple algorithm where,
starting with a perfect crystal structure of $\alpha$-Zr, the Frenkel pairs of
defects are created at random and the microstructure is subsequently relaxed
via direct energy minimisation such that the system evolves purely through the
action of internal stresses arising from the generation of defects. The dose
is measured in units of ‘canonical displacement per atom’ (cdpa) computed as
the ratio of the total number of Frenkel pairs generated by the algorithm to
the number of atoms in the system. CRA simulations assume that vacancies
remain effectively immobile, in turn also immobilising the dislocation part of
the microstructure (Arakawa et al., 2020), leaving the internal fluctuating
stresses as the only remaining factor driving the migration and clustering of
self-interstitial atom defects. In (Warwick et al., 2021) we identified the
approximate temperature and dose rate range where the simulation method
retains its validity when applied to $\alpha$-Zr.

The IIG strains observed in neutron irradiated zirconium alloys with initially
low dislocation densities tend to saturate with increasing dose at
temperatures less than $\sim 300\text{\,}\mathrm{\SIUnitSymbolCelsius}$, see
Adamson et al. (2019). At temperatures above $\sim
300\text{\,}\mathrm{\SIUnitSymbolCelsius}$, saturation persists over shorter
intervals of dose before the strain magnitudes start increasing linearly as a
part of the breakaway growth phenomenon (Holt, 1988). A significant change of
pattern of thermal evolution has also been found above $\sim
300\text{\,}\mathrm{\SIUnitSymbolCelsius}$ in proton irradiated Zircaloy-2
when samples irradiated to $2\text{\,}\mathrm{d}\mathrm{p}\mathrm{a}$ were
annealed for $1\text{\,}\mathrm{h}$ at various temperatures (Topping et al.,
2018). The X-ray line profile measurements performed by Topping et al. (2018)
showed that the a-loop density significantly decreased only at temperatures
above $\sim 300\text{\,}\mathrm{\SIUnitSymbolCelsius}$.

These data offer a valuable insight into the timescales on which thermally
activated processes, including vacancy migration, drive the evolution of
heavily irradiated zirconium. Earlier (Warwick et al., 2021), noting that the
rates of thermally activated processes follow the Arrhenius law (Vineyard,
1957; Landauer and Swanson, 1961; Allnatt and Lidiard, 1993), we showed that
the annealing experiment data imply that the characteristic activation energy
$E_{a}$ for the processes primarily responsible for the observed thermally
activated behaviour must be close to $\sim 2\text{\,}\mathrm{e}\mathrm{V}$.
Also, as described in §2.2, the proton-irradiation defect production dose rate
$\dot{\phi}$ at all depths in our experiments is high and close to $10^{-5}$
dpa s-1. Given this high dose rate, we can estimate an upper bound $\tilde{T}$
on the range of temperatures where the rate of migration of defects stimulated
directly by irradiation is higher than the rate of thermally activated
migration of defects. For a given activation energy $E_{a}$, using the dose
rate model by Nordlund et al. (2018), we find that the two rates are
comparable if

$\displaystyle\dot{\phi}\left({2E_{d}\over
E_{a}}\right)\approx\nu\exp{\left(-\frac{E_{a}}{k_{B}\tilde{T}}\right)},$ (8)

where $E_{d}$ is the threshold displacement energy required for forming a
defect, the attempt frequency is ${\nu\approx\omega_{D}/2\pi}=5.84\times
10^{12}$ s-1, given the Debye frequency ${\omega_{D}=3.67\times
10^{13}\,\mathrm{s}^{-1}}$ (Zarestky, 1979) and $k_{B}=0.861\times 10^{-4}$
eV/K is the Boltzmann constant. Taking $E_{a}=2$ eV and $E_{d}=40$ eV, and
solving equation (8) for $\tilde{T}$, we find $\tilde{T}=625$ K$\approx
350\text{\,}\mathrm{\SIUnitSymbolCelsius}$. Below this temperature, the
eigenrate of thermal relaxation of microstructure is lower than the rate at
which defects are driven by irradiation. Hence, at temperatures below
$\tilde{T}$, the defect structures generated by irradiation evolve
predominantly through fast athermal stress relaxation (Derlet and Dudarev,
2020).

The above esitmate for $\tilde{T}$$\sim
350\text{\,}\mathrm{\SIUnitSymbolCelsius}$ is close to the temperature at
which Topping et al. (2018) observed the occurrence of a significant change in
the thermal response of microstructure during annealing. Notably, the change
of pattern of breakaway growth at $\sim
300\text{\,}\mathrm{\SIUnitSymbolCelsius}$ noted by Holt (1988) was observed
at significantly lower dose rates than those characterising our proton
irradiation experiments. Hence, the above temperature must reflect the
fundamental scale of activation energies associated with microstructural
evolution of Zircaloys under irradiation.

The migration energy of individual vacancies in pure elemental $\alpha$-Zr
(Varvenne et al., 2014) of $\sim 0.5\text{\,}\mathrm{eV}$ is too low to
account for the observed behaviour, and while the thermal diffusion of
vacancies and other point defects affect microstructural evolution, reducing
the overall concentration of defects noted in our analysis, experimental
observations indicate the presence of a rate-limiting process with a higher
activation energy that stabilises the observed dense defect microstructures at
temperatures as high as $300\text{\,}\mathrm{\SIUnitSymbolCelsius}$.

High activation energies are known to be associated with the formation of
immobile vacancy-impurity clusters involving carbon or nitrogen (Fu et al.,
2008; Terentyev et al., 2014; Theodorou et al., 2022). In bcc iron, the
effective migration energy of vacancies is defined by the energy of
dissociation of a cluster involving a vacancy and a carbon dimer (Paxton,
2014), and this dissociation energy can be as high as 2.22 eV (Kabir et al.,
2010), far higher than the activation energy of 0.55 eV characterising vacancy
migration in pure elemental Fe (Fu et al., 2005). Given that the
characteristic formation and migration energies of defects in Zr and Fe are
nearly the same (Dudarev, 2013), the high effective activation energy of the
order of 2 eV seen in experiments on Zr likely result from the impurity effect
similar to that found in Fe. As noted by Kabir et al. (2010), at relatively
low temperatures the vacancy-carbon dimer complexes are immobile, making the
dissociation temperature of these complexes one of the key parameters
determining the response of a material to radiation exposure.

CRA simulations (Derlet and Dudarev, 2020) or the simulations involving the
production of defects by successive collision cascade events (Mason et al.,
2021; Granberg et al., 2023; Boleininger et al., 2023) do not imply the
absence of mobility of defects. Self-interstitial atom defects exhibit the
non-Arrhenius mobility (Dudarev, 2008) and their motion is strongly affected
by elastic strain fields (Dudarev et al., 2010), resulting in the rapid
clustering of these defects into interstitial dislocation loops and,
subsequently, into a dense entangled network of dislocations (Derlet and
Dudarev, 2020; Boleininger et al., 2022). The latter forms spontaneously at
doses above approximately $0.3\text{\,}\mathrm{d}\mathrm{p}\mathrm{a}$ (Mason
et al., 2020). The clustering of self-interstitial defects into dislocation
loops and dislocations stems from the fact that this is a highly energetically
favourable process, releasing up to $E^{f}_{SIA}\approx 3$ eV per self-
interstitial coalescence event (Domain and Legris, 2005; Dudarev, 2013).

The fact that it is the SIA formation energy, fundamentally related to the
strong elastic interaction between the self-interstitial defects, that drives
the evolution of microstructure at relatively low tempeatures rather than the
diffusion of self-interstitial per se, is confirmed by finite-temperature
simulations by Chartier and Marinica (2019). The simulations were performed at
300K and hence included the thermal diffusion of self-interstitial atom
defects, but still exhibited the same pattern of evolution as that predicted
by the CRA simulations (Derlet and Dudarev, 2020). This is confirmed by
experimental observations by Wang et al. (2023) showing the trends similar to
those found in simulations, even though in tungsten the diffusion of self-
interstitial defects occurs at temperatures as low as 27 K (Ehrhart et al.,
1991; Ma and Dudarev, 2019).

The above analysis of ab initio data and experimental information shows that
the temperature interval over which the dynamics of microstructural evolution
changes from the low-temperature mode dominated by microscopic stress
fluctuations (Derlet and Dudarev, 2020) to the high-temperature mode dominated
by the Arrhenius thermally activated diffusion (Allnatt and Lidiard, 1993), in
zirconium alloys spans approximately from
$300\text{\,}\mathrm{\SIUnitSymbolCelsius}$ to
$400\text{\,}\mathrm{\SIUnitSymbolCelsius}$, as illustrated particularly well
by Fig. 8 from Topping et al. (2018). The qualitative picture of
microstructural evolution of zirconium irradiated by high-energy protons can
now be summarised as follows. Proton irradiation produces relatively low
energy recoils, generating defects in the form of Frenkel pairs or small
defect clusters (Boleininger et al., 2023). Self-interstitial atom defects
coalesce into dislocation loops and a dislocation network, whereas vacancies
either diffuse and recombine with the interstitial dislocation loops or
extended dislocation structures, or form immobile vacancy-impurity clusters
(Kabir et al., 2010). These vacancy-impurity clusters immobilise and stabilise
the dislocation microstructure (Arakawa et al., 2020), but dissociate in the
temperature interval from $300\text{\,}\mathrm{\SIUnitSymbolCelsius}$ to
$400\text{\,}\mathrm{\SIUnitSymbolCelsius}$. Over this temperature interval,
the mode of microstructural evolution changes from that dominated by stress
fluctuations and coalescence of self-interstitial defects to the mode
dominated by vacancy diffusion. Over the same interval of temperatures, the
IIG changes from a mode exhibiting saturation to that of runaway growth. Our
observations exhibit the formation of a dense dislocation network, suggesting
that the experimental conditions correspond to the low-temperature rather than
the high-temperature mode of microstructural evolution. The selection of a
simulation approach below reflects and recognises this fact. An algorithm for
modelling microstructural evolution at higher temperatures has to include the
treatment of microscopic stress fluctuations as well as diffusion and
interaction of vacancies and impurities. The development of such an algorithm
remains a challenge for future studies.

The CRA was implemented in the molecular dynamics program lammps (Plimpton,
1995) 111https://lammps.sandia.gov, 3 Mar and 29 Oct 2020 stable builds. For
this study, unless stated otherwise, we present results averaged across all
three Embedded Atom Method (EAM) potentials developed in Ref. (Mendelev and
Ackland, 2007), as is the case in Figure 1. Whilst there are variations
between the potentials with respect to their predicted formation energies and
elastic properties of self-interstitials and vacancies (Mendelev and Ackland,
2007; Varvenne et al., 2014; Varvenne and Clouet, 2017), we have found that
all three potentials qualitatively produce the same macroscopic dimensional
changes and microstructural evolution under the CRA (Warwick et al., 2021).
Furthermore, a similar study also employed the CRA on Zr and predicted the
same trends with a different potential (Tian et al., 2021).

Our simulations employ periodic boundary conditions and supercells containing
$\sim$ 2M and $\sim$ 10M atoms, with the cell edges parallel to the
$[2\bar{1}\bar{1}0]$, $[\bar{1}2\bar{1}0]$ and $[0001]$ directions. Energy
minimisation was performed using a combination of the conjugate gradient and
FIRE algorithms (Bitzek et al., 2006) such that the relaxed force on any atom
was smaller than $1\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\text{\AA}^{-1}$. In
the interest of computational efficiency, the simulation cell shape and size
was kept fixed during relaxation. Whilst these boundary conditions result in a
macroscopic internal stress of $\sim
1\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$, the dimensional changes that would
occur if the cell shape relaxed may nevertheless be accurately computed using
linear elasticity theory; furthermore, the microstructure resulting from
relaxing under zero pressure is similar to that under zero strain (Warwick et
al., 2021; Tian et al., 2021). Dislocations were identified directly from
atomic positions using the Dislocation eXtraction Algorithm (DXA) (Stukowski
et al., 2012). This is achieved by assigning crystal structure types to each
atom using common neighbour analysis (Faken and Jónsson, 1994). Given the
crystal structure of the hcp reference crystal, Burgers circuits are drawn
around regions containing atoms assigned to non-hcp crystal structure in order
to compute Burgers vectors and dislocation lines. Dislocation densities are
computed according to Equation 1.

In order to enable a closer comparison between our simulations and line
profile analysis experiments, we simulated the intensity profile of powder
diffraction patterns for all of the CRA microstructures using the Debye
equation (Debye, 1915) where the intensity of scattered X-rays is proportional
to

$\displaystyle I(q)=\sum_{i,j}{\mathrm{sinc}\left(2\pi qr_{ij}\right)},$ (9)

for wavevector magnitude $q$, and the sum runs over all the pairs of atoms
positioned at $\mathbf{r}_{i}$ and $\mathbf{r}_{j}$ separated by distance
${r_{ij}=|\mathbf{r}_{i}-\mathbf{r}_{j}|}$. Furthermore, in Equation 9 it is
assumed that the atomic form factor is the same for all atoms in the system.
The time required for computing all the pairwise distances $r_{ij}$ for a
system of $N$ atoms scales unfavourably as $N^{2}$ and thus we parallelised
the task. The line profile was calculated without using periodic boundary
conditions and thus the powder was treated as if it were composed of randomly
oriented nano-grains as large as the simulation box. $I(q)$ was computed over
the domain $3\text{\,}\mathrm{n}\mathrm{m}^{-1}$ $\leq q\leq$
$13\text{\,}\mathrm{n}\mathrm{m}^{-1}$ spanning all peaks up to the
$\\{22\bar{4}0\\}$ reflections and the wavenumbers were sampled every
$2\text{\times}{10}^{-3}\text{\,}\mathrm{n}\mathrm{m}^{-1}$.

Data were visualised and processed with the ovito (Stukowski, 2010) and
paraview (Ahrens et al., 2005) software packages. For more details, please
refer to our recent publication (Warwick et al., 2021).

##  3 Results and Discussion

###  3.1 Dislocation structure and distribution

(a)

(b)

(c)

Figure 3: Dislocation microstructures calculated by the CRA simulations using
the MA1 potential with 10M atoms. Green, orange and red dislocations indicate
$\langle 2\bar{1}\bar{1}0\rangle,\langle 1\bar{1}00\rangle$ and unusual
Burgers vectors, respectively. Dots are self-interstitial atoms identified
with the Wigner-Seitz analysis and boundaries of interstitial clusters are
rendered as transparent surfaces. 3(a): Self interstitial point defects and
small dislocation loops at very low dose 0.02 cdpa. 3(b) : Dense dislocation
network at low dose 0.13 cdpa. 3(c) : Extended dislocation network at high
dose 2.00 cdpa with percolating interstitial cluster.

The peak and saturation of dislocation density shown in Figure 1 can be
readily understood by inspecting the spatial configuration of dislocations
generated by our simulations. The nature of defects evolving through internal
stresses in the CRA results in dislocations being almost exclusively formed by
the agglomeration of self-interstitial atoms whilst vacancies remain immobile
and generally form small clusters containing $\mathcal{O}(10)$ vacancies that
are not large enough to relax into dislocation loops. We do find a small
number of much larger clusters containing $\mathcal{O}(100)$ vacancies.
However, these are fundamentally interstitial in origin as the large vacancy
clusters are nothing but vacant spaces in the crystallographic planes formed
by the self-interstitial defects, see Mason et al. (2020); Boleininger et al.
(2023). The renders shown in Figure 3 suggest that the dislocation structure
evolves in three stages. At low dose (Figure 3(a)), small loops form before
coalescing into a dislocation network (Figure 3(b)) at which point the
dislocation density saturates. At high doses (Figure 3(c)), full interstitial-
type atomic planes form and the dislocation network fragments into loops,
resulting in a drop in the dislocation density. The visualised microstructures
were rendered from simulations employing the MA1 potential. When comparing the
interatomic potentials, we discovered that the MA2 and MA3 potentials produce
large populations of twinned regions (Warwick et al., 2021). It seems likely
that these are artefacts of the potentials since such defects are not commonly
observed in experiment. As was noted by Warwick et al. (2021), for the MA3
potential in particular, a large proportion of dislocations coalesce into
these twinned regions whose volume fraction also features a transient peak
followed by saturation. The twinned regions are composed of dense arrays of
dislocations and thus the pattern of evolution of dislocation structures and
saturation in density is common to all three potentials. Thus, the
microstructures derived from MA1 simulations are presented in order to avoid
needlessly complicating our discussion.

Figure 4: Simulated $\\{1\bar{1}02\\}$ X-ray diffraction peak profile as a
function of dose computed from the MA1 microstructures. The variation in the
intensity of the peak and its centre shift is highlighted by the red path.
Application of the cmwp software to our simulated line profiles also shows
that the dislocation density saturates as a function of dose, approaching a
limiting value of $\sim 10^{15}$ m-2.

CRA simulations tend to overestimate the defect content in the high dose limit
(Boleininger et al., 2023), however the qualitative trends observed under a
variety of conditions are predicted accurately (Mason et al., 2020, 2021;
Warwick et al., 2021). This overestimation arises from the lack of re-
crystallisation induced by collision cascades. Producing the primary knock on
atoms with kinetic energies close to the threshold displacement energy results
in defect densities very close to those predicted by the CRA whilst much
larger recoil energies result in defect densities lower by approximately a
factor of approximately 10, see Boleininger et al. (2023). Thus, in agreement
with the analysis by Mason et al. (2020, 2021); Boleininger et al. (2023)
where the defect content was independently assessed using the experimentally
measured deuterium retention, we have scaled down the dislocation density in
Figure 1 by a factor of 10 to enable a direct comparison of our experimental
and simulation data. After applying this scaling it can be seen that the CRA
indeed predicts a qualitatively accurate variation dislocation density
profile. The occurrence of a transient peak of dislocation density at a
moderate dose was also noted earlier in simulations of iron and tungsten
(Chartier and Marinica, 2019; Derlet and Dudarev, 2020).

In order to enable a closer comparison with experiment we simulated a powder
diffraction line profile for all of our CRA microstructures as described in
§2.3. Figure 4 illustrates the evolution of the $\\{1\bar{1}02\\}$ diffraction
peak profile as a function of dose, reflecting the eventual saturation of the
microstructure in the peak profile seen in the limit of high dose. We observe
that in the transient regime, the peak intensity drops before rising and
settling at a steady value. Furthermore, the peak broadens at high dose and
shifts to higher scattering angles, indicating lattice compression. The
formation of extra atomic planes due to the coalescence of dislocation loops
does not result in the volumetric expansion but instead causes lattice
compression because the simulations are performed at constant cell shape and
size. Under zero pressure boundary conditions, we expect the peak centre to
shift to lower wave-numbers instead.

Saturation of the peak profile has been quantified by extracting the
dislocation density from the data using the cmwp software. cmwp is known to
infer dislocation densities that are larger than those determined from TEM
images and this has been attributed to the ability of X-ray line profile
analysis to resolve small loops in power law distributed dislocation loops
(Ungár et al., 2021, 2021). Interestingly, we find that the dislocation
densities computed by cmwp and DXA nevertheless differ by approximately an
order of magnitude although the character of variation of the observed and
simulated dislocation densities as functions of dose are the same. We note
that the difference between the cmwp software and DXA results brings attention
to an important question: at what size is a dislocation loop too small to be
counted as a dislocation object?

Usually, dislocations are considered to be the sources of long-range strain
fluctuations responsible for the X-ray peak line profile broadening detected
in irradiated materials (Wilkens, 1970). On the other hand, small dislocation
loops produce shorter-range strains and in the far-field limit they are
equivalent to point defects, with the associated strain fields resulting in
the Huang diffuse scattering, producing a relatively uniform increase in the
scattered intensity in X-ray diffraction patterns (Simmons and Baluffi, 1958).
Furthermore, when a dislocation loop is so small that the loop diameter is
comparable to the core width of a dislocation then the loop is mostly
comprised of core atoms and its structure can no longer be described using
conventional elasticity (Boleininger et al., 2018; Boleininger and Dudarev,
2019). Thus, one would expect the resulting strains to be unlike those
associated with linear elastic fields of dislocation loops and the core
effects to be significant. Modelling the strain broadening effects associated
with small dislocation loops requires further analysis and we defer it to a
future study. A large proportion of the dislocation objects found in the
simulated microstructures are so small that they should be treated as point
defects by CMWP and would also not be easily detected in TEM images (Zhou et
al., 2006). Determining the nature of such defects could be highly relevant to
determining mechanisms that cause the complex high dose phenomena such as
breakaway growth.

(a)

(b)

Figure 5:  5(a) : Histogram of self-interstitial cluster population sizes
detected in MA1 microstructures containing 10M atoms corresponding to Figures
3. The population sizes of clusters, $N$, follow a power law distribution.
$N_{\mathrm{SIA}}$ denotes the total number of self-interstitial atoms
detected by the Wigner-Seitz algorithm. The data are presented on a log-log
scale. Lines of best fit were calculated using the maximum likelihood
estimates (MLE) for the exponent in a power law distribution stated in Eq. 10.
5(b): Evolution of the MLEs for the power law exponent $\alpha$ as a function
of dose averaged over all three interatomic MA potentials. Shaded regions
represent confidence intervals of the MLE. Data for the doses lower than
$1\text{\,}\mathrm{c}\mathrm{d}\mathrm{p}\mathrm{a}$ are presented on a log-
linear scale whereas in the saturation regime
$>1\text{\,}\mathrm{c}\mathrm{d}\mathrm{p}\mathrm{a}$ a linear-linear scale is
employed. $\alpha$ was averaged over doses larger than
$1\text{\,}\mathrm{c}\mathrm{d}\mathrm{p}\mathrm{a}$ to arrive at the value
$\langle\alpha\rangle=2.28\pm 0.01$ indicated by the horizontal dashed blue
line for a 2M cell, and $\langle\alpha\rangle=2.23\pm 0.01$ (dashed green
line) for a 10M cell.

Size distributions of dislocation loops are often investigated in experiment
(Yi et al., 2015) and so for the purposes of comparison and characterising the
spatial arrangement of dislocations we examined the distribution of defect
cluster - dislocation loop sizes. As the dislocations in this simulation are
mostly of interstitial type, we can gain insight into the statistics of
dislocation structures by examining the statistics of interstitial defect
clusters. Figure 5(a) shows the frequency distribution of clusters containing
$N$ interstitials as calculated by ovito for the doses corresponding to the
renders in Figure 3. For clusters with population sizes $N<10$, the bin widths
of the histogram are equal to 1 whereas at larger cluster population sizes the
bin widths increase logarithmically (Milojević, 2010). The raw data for the
cluster population size frequencies were averaged over the bin widths to
produce the step chart shown. Visually, the histograms appear to follow a
straight line on a log-log scale, suggesting a power law distribution.
Furthermore, simulations employing the CRA have shown evidence for self-
organised critical behaviour (Derlet and Dudarev, 2020) and thus it is
reasonable to test the cluster population size power law distribution
hypothesis such that the probability mass function for clusters containing $N$
interstitials is given by

$p(N)=\frac{1}{\zeta(\alpha)N^{\alpha}}$ (10)

where the Riemann zeta function (Heynsworth and Goldberg, 1965) is defined as
$\zeta(\alpha)=\sum_{k=1}^{\infty}1/k^{\alpha}$ for exponent $\alpha>1$.

We may define an exponent that best fits the data shown in Figure 5(a) via a
maximum likelihood estimation. The likelihood function

$\displaystyle\mathcal{L}\left(\\{N_{i}|i\in[1..N_{\mathrm{tot}}]\\}|p(N|\alpha)\right)=\prod_{i=1}^{N_{\mathrm{tot}}}p(N_{i}|\alpha)$
(11)

returns the probability of observing the $N_{\mathrm{tot}}$ measured data
points $\\{N_{i}|i\in[1..N_{\mathrm{tot}}]\\}$ if they were produced from a
distribution $p(N|\alpha)$ of given parameter(s) $\alpha$. The MLE is produced
by determining the value of $\alpha_{\mathrm{MLE}}$ maximising
$\ln\mathcal{L}$ such that:

$\displaystyle\left.\frac{\partial\ln\mathcal{L}}{\partial\alpha}\right\rvert_{\alpha=\alpha_{\mathrm{MLE}}}=0.$
(12)

Equation 12 was solved numerically using the python package powerlaw (Alstott
et al., 2014) and as shown in Figure 5(b), the MLEs and associated standard
error for $\alpha$ exhibit transients over a range of doses corresponding to
the formation of a dislocation network by the coalescence of loops. The
minimum of this transient appears to be correlated with the peak in
dislocation density before saturating to a constant value at high doses. The
twinned regions that emerge when employing the MA2 and MA3 potentials, as
described at the beginning of §3.1, are also interstitial defect clusters.
Therefore they are included in this statistical analysis and we find the same
trend and remarkably similar values of $\alpha$ across all three MA
potentials. Averaging across the potentials for doses larger then
$1\text{\,}\mathrm{c}\mathrm{d}\mathrm{p}\mathrm{a}$ in the saturation regime,
we find the MLE for exponent $\alpha=2.28\pm 0.01$ for a 2M simulation cell
and $\alpha=2.23\pm 0.01$ for a 10M simulation cell. These values are close
but slightly higher than the exponents found in simulations and observations
of collision cascades in tungsten in the limit of low dose (Sand et al., 2013;
Yi et al., 2015).

Figure 6: Estimation of dislocation density from interstitial cluster
perimeter density. $N$ corresponds to the number of atoms in the cluster and
the perimeter of a given cluster is estimated according to Equation 13.

An immediate observation and the consequence of the power law statistics is
that the overwhelming majority of clusters are small. Thus, even the order of
magnitude of the calculated dislocation density is sensitive to the choice one
makes for a minimum detectable size of a dislocation loop. To show this,
assume that clusters with cluster population sizes $N$ lying in a range
${N_{\text{min}}<N<N_{\text{max}}}$ correspond to a circular platelet of
interstitial atoms i.e. a dislocation loop (Gilbert et al., 2008). The
$N_{\text{min}}$ and $N_{\text{max}}$ thresholds correspond to cluster
population sizes that are either sufficiently small to be treated as point
defects, or are large enough to be close to the threshold for forming a
percolating atomic plane consisting of interconnected dislocation loops.
Assume that each platelet contains one extra atom per atomic string in the
direction of the Burgers vector (Boleininger et al., 2018; Boleininger and
Dudarev, 2019), and involves $N$ atoms with volume $\Omega$ equal to the
atomic volume of $\alpha$-Zr. Using the formula $({\bf b}\cdot{\bf A})$ for
the volume of a dislocation loop, its perimeter can then be estimated as

$P_{S}=2\pi\sqrt{{N\Omega\over\pi b}},$ (13)

where $b$ is the length of the Burgers vector along the normal to the loop
habit plane. Here, we assume that the platelets form $a$-type dislocation
loops such that $b$ is equal to the $a$ lattice parameter.

Summing the perimeters corresponding to cluster population size distributions
between various choices of $N_{\text{min}}$ and $N_{\text{max}}$ provides a
measure of the difference in dislocation density due to different choices of
thresholds. Typically, a dislocation core extends over about five interatomic
distances (Boleininger et al., 2018; Boleininger and Dudarev, 2019) and thus
one would reasonably argue that $N_{\text{min}}$ should at least be greater
than five. Furthermore, we should not include the fully formed planes
corresponding to a cluster percolating through the periodic boundaries, thus
setting $N_{\text{max}}=10^{3}$. In Figure 6 we observe an order of magnitude
difference between the results involving the counting of all the viable
clusters smaller than the population size of the percolating cluster
($5<N<10^{3}$), and the values found by counting all the clusters but
increasing the lower threshold to exclude clusters containing less than 100
self-interstitial atoms. This example illustrates how the power law
distribution of defect sizes in highly irradiated microstructures can
profoundly affect how dislocations are counted and their density quantified.

We may also employ Equation 13 to estimate the power law exponent for the
distribution of loops with respect to loop diameter $D$ where it is apparent
that, assuming that all the loops have the same Burgers vector $b$, $N\propto
D^{2}$. Treating $N$ as a continuous variable, we derive the probability
density function as a function of loop diameter

$\displaystyle p_{D}(D)=\left\lvert\frac{\mathrm{d}N}{\mathrm{d}D}\right\rvert
p(N(D))\propto\frac{1}{D^{2\alpha-1}},$ (14)

where $\alpha$ is the exponent entering Equation 10. Hence the diameter of
circular loops is expected to be power law distributed with exponent
$\beta=2\alpha-1$. Using the data derived from the CRA simulations in the high
dose limit illustrated in Figure 5(b) for a 10M atom cell, we find that
$\beta=3.46$. In experiments performed by Ungár et al. (2021) dislocation
diameter data measured by TEM were combined with an XRD line profile
measurement of the total dislocation density. When fitting the dislocation
size distribution to a power law the exponent was found to lie in the range
$3\leq\beta\leq 4$, which agrees well with the values derived from our
simulations.

Concluding this section, we note that power law exponent values $\alpha=2$ and
$\beta=3$ appear to represent natural low limits characterising the power law
statistics of populations of dislocation loops in a heavily irradiated
material. Indeed, any power law distribution with $\alpha<2$ would imply a
divergent total count of interstitial defects contained in the loops, a
paradox that can only be resolved by recognising that loops of very large size
are nothing but extra crystallographic atomic planes, which through elastic
interactions would tend to modify the finite-size part of the dislocation loop
distribution so as to make it normalisable, in this way steering the value of
exponent $\alpha$ towards and above the limiting value of 2.

###  3.2 Stored energy

Whilst analysing the configuration of atoms helps to describe how the
microstructure evolves, on its own this approach provides limited insight into
why the observed processes take place. Evidently, the repeated creation and
relaxation of defects forces the microstructure into a state that has the
energy higher than that of a single crystal. By examining the excess energy
$E_{exc}$ that the system is driven to, we may determine how the accumulation
of point defects produces the observed dislocation density. For an irradiated
microstructure at dose $\phi$ containing $N$ atoms with total energy
$E_{total}(\phi)$ we may measure the excess (stored) energy as

$\displaystyle E_{exc}(\phi):=E_{total}(\phi)-NE_{coh},$ (15)

where $E_{coh}$ is the cohesive energy of $\alpha$-Zr.

Our simulations are carried out under zero global strain boundary conditions.
Whilst this is computationally convenient, in reality specimens are often
irradiated under zero applied stress, allowing the body as a whole to undergo
strain. Assuming linear elasticity, we may correct for this and remove the
stored elastic energy $E_{el}$ from our simulation results. The defects
produced by radiation damage induce eigenstrains, also known as residual
strains, that act as sources of elastic strain. Let $\epsilon_{ij}^{0}$ denote
the elastic strain that would come about under zero applied stress boundary
conditions. In this case, the potential energy of the body is lowered by doing
work equal to $\frac{1}{2}\sigma_{ij}\epsilon_{ij}^{0}$ where Hooke’s law
determines the stress $\sigma_{ij}=C_{ijkl}\epsilon_{kl}^{0}$ and the elastic
constants tensor is $C_{ijkl}$. Enforcing zero global strain requires that the
integral of the strain $\epsilon_{ij}$ over the body with volume $V$ is zero
or, equivalently, that one has zero volume averaged strain

$\displaystyle\langle\epsilon_{ij}\rangle:=\frac{1}{V}\int_{V}\mathrm{d}^{3}\mathbf{x}^{\prime}\,\epsilon_{ij}(\mathbf{x}^{\prime})=0.$
(16)

This boundary condition is satisfied by the strain

$\displaystyle\epsilon_{ij}(\mathbf{x})=\epsilon_{ij}^{0}(\mathbf{x})-\langle\epsilon^{0}_{ij}\rangle.$
(17)

Applying Hooke’s law to Equation 17, we observe that the system is under a
state of global stress
$\langle\sigma_{ij}\rangle=C_{ijkl}\langle\epsilon^{0}_{kl}\rangle$ and we may
treat this as the stress that develops in our simulations. Therefore, we can
compute the stored elastic energy

$\displaystyle E_{el}=\frac{1}{2}\langle\sigma_{ij}\rangle
S_{ijkl}\langle\sigma_{kl}\rangle,$ (18)

where $S_{ijkl}$ is the elastic compliance tensor (Warwick et al., 2021)
related to $C_{ijkl}$ by
$S_{ijpq}C_{pqkl}=\frac{1}{2}\left(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right)$.
We find across all the employed potentials that $E_{el}$ accounts only for a
small fraction $<10\%$ of the stored energy, meaning that the vast majority of
stored energy is contained in the point defect centres, dislocation cores and
fluctuations arising from the elastic fields of these defects. Indeed, when a
high dose snapshot was explicitly relaxed under zero pressure we found that
the difference in excess energy to that corrected for by our estimate of
$E_{el}$ is of the same order of magnitude. Furthermore, the microstructure
did not change significantly providing more indication of the minor role of
the global elastic energy.

(a)

(b)

(c)

Figure 7: Contributions of microstructural defects to the excess energy.
Results shown are for the MA1 potential using 10 million atoms. 7(a):
Comparison of vacancy energy $E_{vac}$ estimated by assuming all $N_{V}$
vacancies are isolated with formation energy $E^{f}_{1vac}$ (blue curve) with
$E_{vac}$ calculated via explicitly relaxing pristine $\alpha$-Zr containing
the same number and arrangement of vacancies. 7(b) : Breakdown of total excess
energy $E_{exc}$ into elastic $E_{el}$, interstitial $E_{int}$ and vacancy
$E_{vac}$ contributions. 7(c): Comparison of interstitial excess energy,
concentration and dislocation density that exhibits the similarity in profile
between all three quantities. A peak occurs in all three profiles near
$0.1\text{\,}\mathrm{c}\mathrm{d}\mathrm{p}\mathrm{a}$.

The remaining stored energy is confined in the two classes of defects
associated with vacancies and interstitials. Given that vacancies do not
cluster significantly, we may estimate their energetic contribution as

$\displaystyle E_{vac}(\phi)=N_{vac}(\phi)E^{f}_{1vac},$ (19)

where, for dose $\phi$, the number of vacancies is denoted $N_{vac}$ and
$E^{f}_{1vac}$ is the formation energy of a single vacancy. Values of
$E^{f}_{1vac}$ for each of the potentials used in this study may be found in
Mendelev and Ackland (2007). To check the validity of Equation 19 we isolated
the $N_{vac}$ vacancies identified by Wigner-Seitz analysis and subsequently
relaxed an $N$-atom supercell of pristine $\alpha$-Zr containing the same
number of vacancies in the same positions. The simulation cell was relaxed and
the formation energy was computed as

$\displaystyle E_{vac}=E^{total}_{vac}-(N-N_{vac})E_{coh},$ (20)

where $E^{total}_{vac}$ is the resulting total energy. Figure 7(a) summarises
this process and shows the resulting formation energy for these arrangements
of vacancies. Thus we observe that Equation 19 is a good approximation,
providing further evidence that the majority of $E_{vac}$ is due to isolated
vacancies.

The contribution to $E_{exc}$ from the formation energy of small interstitial
clusters and dislocation content may now be calculated as

$\displaystyle E_{int}=E_{exc}-E_{el}-E_{vac}.$ (21)

In Figure 7(b) we show the relative proportion of elastic, interstitial and
vacancy contributions to the total excess energy where it is evident that the
elastic contribution is in the minority, vacancy contributions dominate and
the interstitial contribution follows a similar trend to the dislocation
density profile shown in Figure 1. The profile of $E_{int}$ is shown in Figure
7(c) together with the interstitial concentration and dislocation density in
order to highlight their correlation with each other.

Recently, Differential Scanning Calorimetry (DSC) experiments were performed
to measure the stored energy of irradiated titanium as a means of inferring
the number of defects present (Hirst et al., 2022). The measurements indicated
stored energies associated with irradiation induced defects to be on the order
of $0.1\text{\,}\mathrm{J}\mathrm{g}^{-1}$. Their analysis also allowed the
authors to infer the presence of defects that are invisible to TEM imaging. At
$0.07\text{\,}\mathrm{c}\mathrm{d}\mathrm{p}\mathrm{a}$ in Figure 7(c), we
find that the peak in specific energy associated with $E_{int}$ is
$25\text{\,}\mathrm{J}\mathrm{g}^{-1}$ that subsequently drops and plateaus at
high dose to $12\text{\,}\mathrm{J}\mathrm{g}^{-1}$. At elevated temperatures,
the defect content is likely to be $\sim 10\%$ of that calculated in our
simulations cf. (Mason et al., 2020, 2021) and thus we expect the associated
stored energy in Zr to be comparable to $1\text{\,}\mathrm{J}\mathrm{g}^{-1}$.

##  4 Conclusions

In summary, we have performed experiments and simulations showing that the
dislocation density in irradiated zirconium and zircaloys exhibits a peak at a
moderate dose and then saturates at doses greater than
$1\text{\,}\mathrm{d}\mathrm{p}\mathrm{a}$. Simulations indicate that this
occurs in a regime of dose rate and temperature where microstructural
evolution is predominantly driven by stress relaxation. The material enters a
critical state at $\sim 1\text{\,}\mathrm{d}\mathrm{p}\mathrm{a}$, where
interstitial clusters grow to a sufficient size to percolate the volume of the
material. At high dose, the population of smaller clusters and dislocation
loops is distributed as a function of cluster defect content according to a
power law statistics with the exponent close to $\alpha\approx 2.2$. As a
function of defect diameter, this results in a power law distribution of
defect clusters with exponent of $\beta\approx 3.5$, which compares favourably
with the range of values $3\leq\beta\leq 4$ derived from experimental
observations (Ungár et al., 2021).The analysis highlights the significance of
precise definition of defect sizes included in the measured dislocation
densities. Irrespectively of the statistics of dislocation structures, the
trend in the dislocation density evolution in zirconium irradiated at
temperatures below $\sim 350^{\circ}\>\mathrm{C}$ is clear, the dislocation
density saturates as a function of dose.

##  5 Acknowledgements

This work received funding from the RCUK Energy Programme Grant No.
EP/W006839/1 and MIDAS EPSRC Grant No. EP/S01702X/1, and was partially 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
authors 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. We acknowledge DESY (Hamburg, Germany), a
member of the Helmholtz Association HGF, for the provision of experimental
facilities. We gratefully acknowledge the use of the high-performance
computing facility MARCONI (Bologna, Italy) provided by EUROfusion, and
computing resources supplied by the IRIS (STFC) Consortium. This work also
received support from the EPSRC Access to HPC Programme on the ARCHER2 UK
National Supercomputing Service (http://www.archer2.ac.uk).

## References

  * Adamson et al. (2019) Adamson, R.B., Coleman, C.E., Griffiths, M., 2019.  Irradiation creep and growth of zirconium alloys: A critical review.  Journal of Nuclear Materials 521, 167–244.  doi:10.1016/j.jnucmat.2019.04.021.
  * Ahrens et al. (2005) Ahrens, J., Geveci, B., Law, C., 2005.  ParaView: An End-User Tool for Large-Data Visualization, in: Hansen, C.D., Johnson, C.R. (Eds.), Visualization Handbook. Elsevier, pp. 717–731.  doi:10.1016/B978-012387582-2/50038-1.
  * Allnatt and Lidiard (1993) Allnatt, A.R., Lidiard, A.B., 1993\.  Atomic Transport in Solids.  Cambridge University Press, Cambridge, England.  doi:10.1017/CBO9780511563904.
  * Alstott et al. (2014) Alstott, J., Bullmore, E., Plenz, D., 2014.  powerlaw: A Python Package for Analysis of Heavy-Tailed Distributions.  PLoS ONE 9, e85777.  doi:10.1371/journal.pone.0085777.
  * Arakawa et al. (2020) Arakawa, K., Marinica, M.C., Fitzgerald, S., Proville, L., Nguyen-Manh, D., Dudarev, S.L., Ma, P.W., Swinburne, T.D., Goryaeva, A.M., Yamada, T., Amino, T., Arai, S., Yamamoto, Y., Higuchi, K., Tanaka, N., Yasuda, H., Yasuda, T., Mori, H., 2020\.  Quantum de-trapping and transport of heavy defects in tungsten.  Nature Materials 19, 508\.  doi:10.1038/s41563-019-0584-0.
  * Arsenlis and Parks (1999) Arsenlis, A., Parks, D., 1999\.  Crystallographic aspects of geometrically-necessary and statistically-stored dislocation density.  Acta Materialia 47, 1597–1611.  doi:10.1016/S1359-6454(99)00020-8.
  * Balogh et al. (2016) Balogh, L., Long, F., Daymond, M.R., 2016.  Contrast factors of irradiation-induced dislocation loops in hexagonal materials.  Journal of Applied Crystallography 49, 2184–2200.  doi:10.1107/S1600576716018136.
  * Bitzek et al. (2006) Bitzek, E., Koskinen, P., Gähler, F., Moseler, M., Gumbsch, P., 2006.  Structural Relaxation Made Simple.  Physical Review Letters 97, 170201\.  doi:10.1103/PhysRevLett.97.170201.
  * Boleininger and Dudarev (2019) Boleininger, M., Dudarev, S.L., 2019\.  Continuum model for the core of a straight mixed dislocation.  Physical Review Materials 3, 093801\.  doi:10.1103/PhysRevMaterials.3.093801.
  * Boleininger et al. (2022) Boleininger, M., Dudarev, S.L., Mason, D.R., Martínez, E., 2022\.  Volume of a dislocation network.  Physical Review Materials 6, 063601\.  doi:10.1103/PhysRevMaterials.6.063601.
  * Boleininger et al. (2023) Boleininger, M., Mason, D.R., Sand, A.E., Dudarev, S.L., 2023\.  Microstructure of a heavily irradiated metal exposed to a spectrum of atomic recoils.  Scientific Reports 13, 1684\.  doi:10.1038/s41598-022-27087-w.
  * Boleininger et al. (2018) Boleininger, M., Swinburne, T.D., Dudarev, S.L., 2018.  Atomistic-to-continuum description of edge dislocation core: Unification of the Peierls-Nabarro model with linear elasticity.  Physical Review Materials 2, 083803\.  doi:10.1103/PhysRevMaterials.2.083803.
  * Chartier and Marinica (2019) Chartier, A., Marinica, M.C., 2019\.  Rearrangement of interstitial defects in alpha-Fe under extreme condition.  Acta Materialia 180, 141–148.  doi:10.1016/j.actamat.2019.09.007.
  * Choi and Kim (2013) Choi, S.I., Kim, J.H., 2013\.  Radiation-induced dislocation and growth behavior of zirconium and zirconium alloys - a review.  Nuclear Engineering and Technology 45, 385–392.  doi:10.5516/NET.07.2013.035.
  * Christensen et al. (2020) Christensen, M., Wolf, W., Freeman, C., Wimmer, E., Adamson, R., Griffiths, M., Mader, E., 2020.  Vacancy loops in Breakaway Irradiation Growth of zirconium: Insight from atomistic simulations.  Journal of Nuclear Materials 529, 151946.  doi:10.1016/j.jnucmat.2019.151946.
  * Das et al. (2018) Das, S., Hofmann, F., Tarleton, E., 2018.  Consistent determination of geometrically necessary dislocation density from simulations and experiments.  International Journal of Plasticity 109, 18–42.  doi:10.1016/j.ijplas.2018.05.001.
  * Debye (1915) Debye, P., 1915.  Zerstreuung von Röntgenstrahlen.  Annalen der Physik 351, 809–823.  doi:10.1002/andp.19153510606.
  * Derlet and Dudarev (2020) Derlet, P.M., Dudarev, S.L., 2020\.  Microscopic structure of a heavily irradiated material.  Physical Review Materials 4, 023605\.  doi:10.1103/PhysRevMaterials.4.023605.
  * Domain and Legris (2005) Domain, C., Legris, A., 2005\.  Ab initio atomic-scale determination of point-defect structure in hcp zirconium.  Philosophical Magazine 85, 569–575.  doi:10.1080/14786430412331334625.
  * Dudarev (2008) Dudarev, S.L., 2008.  The non-Arrhenius migration of interstitial defects in bcc transition metals.  Comptes Rendus Physique 9, 409 – 417.  doi:10.1016/j.crhy.2007.09.019.
  * Dudarev (2013) Dudarev, S.L., 2013.  Density functional theory models for radiation damage.  Annu. Rev. Mater. Res. 43, 35 – 61.  doi:10.1146/annurev-matsci-071312-121626.
  * Dudarev et al. (2010) Dudarev, S.L., Gilbert, M.R., Arakawa, K., Mori, H., Yao, Z., Jenkins, M.L., Derlet, P.M., 2010.  Langevin model for real-time Brownian dynamics of interacting nanodefects in irradiated metals.  Physical Review B 81, 224107\.  doi:10.1103/PhysRevB.81.224107.
  * Dudarev and Ma (2018) Dudarev, S.L., Ma, P.W., 2018\.  Elastic fields, dipole tensors, and interaction between self-interstitial atom defects in bcc transition metals.  Physical Review Materials 2, 033602\.  doi:10.1103/PhysRevMaterials.2.033602.
  * Ehrhart et al. (1991) Ehrhart, P., Jung, P., Schultz, H., Ullmaier, H., 1991\.  Landolt-Börnstein - Group III Condensed Matter · Volume 25: “Atomic Defects in Metals”.  Springer-Verlag Berlin Heidelberg.  doi:10.1007/10011948_45.
  * Faken and Jónsson (1994) Faken, D., Jónsson, H., 1994\.  Systematic analysis of local atomic structure combined with 3D computer graphics.  Computational Materials Science 2, 279–286.  doi:10.1016/0927-0256(94)90109-0.
  * Fu et al. (2005) Fu, C.C., Dalla Torre, J., Willaime, F., Bocquet, J.L., Barbu, A., 2005.  Multiscale modelling of defect kinetics in irradiated iron.  Nature Materials 4, 68 – 74.  doi:10.1038/nmat1286.
  * Fu et al. (2008) Fu, C.C., Meslin, E., Barbu, A., Willaime, F., Oison, V., 2008.  Effect of C on vacancy migration in $\alpha$-iron.  Solid State Phenomena 139, 157 – 164.  doi:10.4028/www.scientific.net/SSP.139.157.
  * Gilbert et al. (2008) Gilbert, M.R., Dudarev, S.L., Derlet, P.M., Pettifor, D.G., 2008\.  Structure and metastability of mesoscopic vacancy and interstitial loop defects in iron and tungsten.  J. Phys.: Condens. Matter 20, 345214\.  doi:10.1088/0953-8984/20/34/345214.
  * Granberg et al. (2023) Granberg, F., Mason, D., Byggmästar, J., 2023.  Effect of simulation technique on the high-dose damage in tungsten.  Computational Materials Science 217, 111902.  doi:10.1016/j.commatsci.2022.111902.
  * Griffiths (2020) Griffiths, M., 2020.  1.11 - Irradiation Growth, in: Konings, R.J., Stoller, R.E. (Eds.), Comprehensive Nuclear Materials. second ed.. Elsevier, Oxford. volume 1, pp. 367–405.  doi:10.1016/B978-0-12-803581-8.11646-7.
  * Groma and Borbély (2004) Groma, I., Borbély, A., 2004\.  X-ray Peak Broadening Due to Inhomogeneous Dislocation Distributions, in: Mittemeijer, E.J., Scardi, P. (Eds.), Diffraction Analysis of the Microstructure of Materials. Springer, Berlin, Heidelberg. Springer Series in Materials Science, pp. 287–307.  doi:10.1007/978-3-662-06723-9_11.
  * Heynsworth and Goldberg (1965) Heynsworth, E.Y., Goldberg, K., 1965\.  Bernoulli and Euler Polynomials, Riemann Zeta Function, in: Abramovitz, M., Stegun, I. (Eds.), Handbook of Mathematical Functions. Dover, New York, pp. 803 – 819.
  * Hirst et al. (2022) Hirst, C.A., Granberg, F., Kombaiah, B., Cao, P., Middlemas, S., Kemp, R.S., Li, J., Nordlund, K., Short, M.P., 2022.  Revealing hidden defects through stored energy measurements of radiation damage.  Science Advances 8.  doi:10.1126/sciadv.abn2733.
  * Holt (1988) Holt, R., 1988.  Mechanisms of irradiation growth of alpha-zirconium alloys.  Journal of Nuclear Materials 159, 310–338.  doi:10.1016/0022-3115(88)90099-2.
  * Hull and Bacon (2011) Hull, D., Bacon, D., 2011.  Introduction to Dislocations : Chapter 1 - Defects in Crystals.  fifth ed., Butterworth-Heinemann, Oxford.  doi:10.1016/B978-0-08-096672-4.00001-3.
  * Jones et al. (2016) Jones, R.E., Zimmerman, J.A., Po, G., 2016.  Comparison of Dislocation Density Tensor Fields Derived from Discrete Dislocation Dynamics and Crystal Plasticity Simulations of Torsion.  Journal of Materials Science Research 5, 44.  doi:10.5539/jmsr.v5n4p44.
  * Kabir et al. (2010) Kabir, M., Lau, T.T., Lin, X., Yip, S., Van Vliet, K.J., 2010\.  Effects of vacancy-solute clusters on diffusivity in metastable Fe-C alloys.  Physical Review B 82, 134112\.  doi:10.1103/PhysRevB.82.134112.
  * Kamminga and Delhez (2000) Kamminga, J.D., Delhez, R., 2000\.  Calculation of diffraction line profiles from specimens with dislocations. A comparison of analytical models with computer simulations.  Journal of Applied Crystallography 33, 1122–1127.  doi:10.1107/S0021889800006750.
  * Landauer and Swanson (1961) Landauer, R., Swanson, J.A., 1961\.  Frequency factors in the thermally activated processes.  Physical Review 121, 1668 – 1674.  doi:10.1103/PhysRev.121.1668.
  * Lemaignan (2012) Lemaignan, C., 2012.  2.07 - Zirconium Alloys: Properties and Characteristics, in: Konings, R.J.M. (Ed.), Comprehensive Nuclear Materials. Elsevier. volume 2, pp. 217–232.  doi:10.1016/B978-0-08-056033-5.00015-X.
  * Ma and Dudarev (2019) Ma, P.W., Dudarev, S.L., 2019\.  Symmetry-broken self-interstitial defects in chromium, molybdenum, and tungsten.  Physical Review Materials 3, 043606\.  doi:10.1103/PhysRevMaterials.3.043606.
  * Mandadapu et al. (2014) Mandadapu, K.K., Jones, R.E., Zimmerman, J.A., 2014.  On the microscopic definitions of the dislocation density tensor.  Mathematics and Mechanics of Solids 19, 744–757.  doi:10.1177/1081286513486792.
  * Mason et al. (2020) Mason, D.R., Das, S., Derlet, P.M., Dudarev, S.L., London, A.J., Yu, H., Phillips, N.W., Yang, D., Mizohata, K., Xu, R., Hofmann, F., 2020.  Observation of Transient and Asymptotic Driven Structural States of Tungsten Exposed to Radiation.  Physical Review Letters 125, 225503\.  doi:10.1103/PhysRevLett.125.225503.
  * Mason et al. (2021) Mason, D.R., Granberg, F., Boleininger, M., Schwarz-Selinger, T., Nordlund, K., Dudarev, S.L., 2021\.  Parameter-free quantitative simulation of high-dose microstructure and hydrogen retention in ion-irradiated tungsten.  Physical Review Materials 5, 095403\.  doi:10.1103/PhysRevMaterials.5.095403.
  * Mendelev and Ackland (2007) Mendelev, M.I., Ackland, G.J., 2007\.  Development of an interatomic potential for the simulation of phase transformations in zirconium.  Philosophical Magazine Letters 87, 349–359.  doi:10.1080/09500830701191393.
  * Milojević (2010) Milojević, S., 2010.  Power law distributions in information science: Making the case for logarithmic binning.  Journal of the American Society for Information Science and Technology 61, 2417–2425.  doi:10.1002/asi.21426.
  * Nicodemus and Staub (1953) Nicodemus, D.B., Staub, H.H., 1953\.  Fission Neutron Spectrum of U235.  Physical Review 89, 1288\.  doi:10.1103/PhysRev.89.1288.
  * Nordlund et al. (2018) Nordlund, K., Zinkle, S.J., Sand, A.E., Granberg, F., Averback, R.S., Stoller, R., Suzudo, T., Malerba, L., Banhart, F., Weber, W.J., Willaime, F., Dudarev, S.L., Simeone, D., 2018.  Improving atomic displacement and replacement calculations with physically realistic damage models.  Nature Communications 9, 1084\.  doi:10.1038/s41467-018-03415-5.
  * Nye (1953) Nye, J., 1953.  Some geometrical relations in dislocated crystals.  Acta Metallurgica 1, 153–162.  doi:10.1016/0001-6160(53)90054-6.
  * Onimus and Bechade (2012) Onimus, F., Bechade, J.L., 2012\.  4.01 - Radiation Effects in Zirconium Alloys, in: Konings, R.J. (Ed.), Comprehensive Nuclear Materials. Elsevier. volume 4, pp. 1–31.  doi:10.1016/B978-0-08-056033-5.00064-1.
  * Onimus et al. (2022) Onimus, F., Gélébart, L., Brenner, R., 2022.  Polycrystalline simulations of in-reactor deformation of recrystallized Zircaloy-4 tubes: Fast Fourier Transform computations and mean-field self-consistent model.  International Journal of Plasticity 153, 103272.  doi:10.1016/j.ijplas.2022.103272.
  * Paxton (2014) Paxton, A.T., 2014.  From quantum mechanics to physical metallurgy of steels.  Materials Science and Technology 30, 1063.  doi:10.1179/1743284714Y.0000000521.
  * Plimpton (1995) Plimpton, S., 1995.  Fast Parallel Algorithms for Short-Range Molecular Dynamics.  Journal of Computational Physics 117, 1–19.  doi:10.1006/jcph.1995.1039.
  * Pomerance (1951) Pomerance, H., 1951.  Thermal Neutron Capture Cross Sections.  Physical Review 83, 641–645.  doi:10.1103/PhysRev.83.641.
  * Ribárik et al. (2020) Ribárik, G., Jóni, B., Ungár, T., 2020.  The Convolutional Multiple Whole Profile (CMWP) Fitting Method, a Global Optimization Procedure for Microstructure Determination.  Crystals 10, 623\.  doi:10.3390/cryst10070623.
  * Rickover et al. (1975) Rickover, H.G., Geiger, L.D., Lustman, B., 1975.  Technical Report. Technical Information Center. U.S. Department of Energy.  doi:10.2172/4240391.
  * Sand et al. (2013) Sand, A.E., Dudarev, S.L., Nordlund, K., 2013.  High-energy collision cascades in tungsten: Dislocation loops structure and clustering scaling laws.  EPL 103, 46003.  doi:10.1209/0295-5075/103/46003.
  * Schoeck (1962) Schoeck, G., 1962.  Correlation between Dislocation Length and Density.  Journal of Applied Physics 33, 1745–1747.  doi:10.1063/1.1728821.
  * Sears (2006) Sears, V.F., 2006.  Neutron scattering lengths and cross sections.  Neutron News 3, 26–37.  doi:10.1080/10448639208218770.
  * Seymour et al. (2017) Seymour, T., Frankel, P., Balogh, L., Ungár, T., Thompson, S., Jädernäs, D., Romero, J., Hallstadius, L., Daymond, M., Ribárik, G., Preuss, M., 2017.  Evolution of dislocation structure in neutron irradiated Zircaloy-2 studied by synchrotron X-ray diffraction peak profile analysis.  Acta Materialia 126, 102–113.  doi:10.1016/j.actamat.2016.12.031.
  * Simmons and Baluffi (1958) Simmons, R.O., Baluffi, R.W., 1958\.  X-Ray Study of Deuteron-Irradiated Copper near 10∘K.  Physical Review 109, 1142\.  doi:10.1103/PhysRev.109.1142.
  * Stoller et al. (2013) Stoller, R., Toloczko, M., Was, G., Certain, A., Dwaraknath, S., Garner, F., 2013\.  On the use of SRIM for computing radiation damage exposure.  Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 310, 75–80.  doi:10.1016/j.nimb.2013.05.008.
  * Stukowski (2010) Stukowski, A., 2010.  Visualization and analysis of atomistic simulation data with OVITO-the Open Visualization Tool.  Modelling and Simulation in Materials Science and Engineering 18, 015012.  doi:10.1088/0965-0393/18/1/015012.
  * Stukowski et al. (2012) Stukowski, A., Bulatov, V.V., Arsenlis, A., 2012.  Automated identification and indexing of dislocations in crystal interfaces.  Modelling and Simulation in Materials Science and Engineering 20, 085007.  doi:10.1088/0965-0393/20/8/085007.
  * Terentyev et al. (2014) Terentyev, D., Heinola, K., Bakaev, A., Zhurkin, E.E., 2014\.  Carbon–vacancy interaction controls lattice damage recovery in iron.  Scripta Materialia 86, 9 – 12.  doi:10.1016/j.scriptamat.2014.04.003.
  * Theodorou et al. (2022) Theodorou, A., Syskaki, M.A., Kotsina, Z., Axiotis, M., Apostolopoulos, G., Fu, C.C., 2022\.  Interactions between irradiation defects and nitrogen in $\alpha$-Fe: an integrated experimental and theoretical study.  Acta Materialia 239, 118227\.  doi:10.1016/j.actamat.2022.118227.
  * Tian et al. (2021) Tian, J., Wang, H., Feng, Q., Zheng, J., Liu, X., Zhou, W., 2021.  Heavy radiation damage in alpha zirconium at cryogenic temperature: A computational study.  Journal of Nuclear Materials 555, 153159.  doi:10.1016/j.jnucmat.2021.153159.
  * Topping et al. (2018) Topping, M., Ungár, T., Race, C.P., Harte, A., Garner, A., Baxter, F., Dumbill, S., Frankel, P., Preuss, M., 2018.  Investigating the thermal stability of irradiation-induced damage in a zirconium alloy with novel in situ techniques.  Acta Materialia 145, 255–263.  doi:10.1016/j.actamat.2017.11.051.
  * Ungár et al. (2021) Ungár, T., Frankel, P., Ribárik, G., Race, C.P., Preuss, M., 2021.  Size-distribution of irradiation-induced dislocation-loops in materials used in the nuclear industry.  Journal of Nuclear Materials 550, 152945.  doi:10.1016/j.jnucmat.2021.152945.
  * Ungár et al. (2021) Ungár, T., Ribarik, G., Topping, M., Jones, R.M.A., Xu, X.D., Hulse, R., Harte, A., Tichy, G., Race, C.P., Frankel, P., Preuss, M., 2021.  Characterizing dislocation loops in irradiated polycrystalline Zr alloys by X-ray line profile analysis of powder diffraction patterns with satellites.  Journal of Applied Crystallography 54, 803–821.  doi:10.1107/S1600576721002673.
  * Ungár et al. (1999) Ungár, T., Dragomir, I., Révész, Á., Borbély, A., 1999\.  The contrast factors of dislocations in cubic crystals: the dislocation model of strain anisotropy in practice.  Journal of Applied Crystallography 32, 992–1002.  doi:10.1107/S0021889899009334.
  * Varvenne and Clouet (2017) Varvenne, C., Clouet, E., 2017\.  Elastic dipoles of point defects from atomistic simulations.  Physical Review B 96, 224103\.  doi:10.1103/PhysRevB.96.224103.
  * Varvenne et al. (2014) Varvenne, C., Mackain, O., Clouet, E., 2014.  Vacancy clustering in zirconium: An atomic-scale study.  Acta Materialia 78, 65–77.  doi:10.1016/j.actamat.2014.06.012.
  * Vineyard (1957) Vineyard, G.H., 1957.  Frequency factors and isotope effects in solid state rate processes.  J. Phys. Chem. Solids 3, 121 – 127.  doi:10.1016/0022-3697(57)90059-8.
  * Wang et al. (2023) Wang, S., Guo, W., Schwarz-Selinger, T., Yuan, Y., Ge, L., Cheng, L., Zhang, X., Cao, X., Fu, E., Lu, G.H., 2023.  Dynamic equilibrium of displacement damage defects in heavy-ion irradiated tungsten.  Acta Materialia , 118578doi:10.1016/j.actamat.2022.118578.
  * Warwick et al. (2021) Warwick, A.R., Boleininger, M., Dudarev, S.L., 2021.  Microstructural complexity and dimensional changes in heavily irradiated zirconium.  Physical Review Materials 5, 113604\.  doi:10.1103/PhysRevMaterials.5.113604.
  * Wilkens (1970) Wilkens, M., 1970.  Theoretical Aspects of Kinematical X-ray Diffraction Profiles from Crystals Containing Dislocation Distribution, in: Simmons, J.A., de Wit, R., Bullough, R. (Eds.), Fundamental Aspects of Dislocation Theory, Washington, D.C.. pp. 1195–1221.
  * Willaime and Massobrio (1989) Willaime, F., Massobrio, C., 1989\.  Temperature-Induced hcp-bcc Phase Transformation in Zirconium: A Lattice and Molecular-Dynamics Study Based on an N-Body Potential.  Physical Review Letters 63, 2244–2247.  doi:10.1103/PhysRevLett.63.2244.
  * Yi et al. (2015) Yi, X., Sand, A.E., Mason, D.R., Kirk, M.A., Roberts, S.G., Nordlund, K., Dudarev, S.L., 2015\.  Direct observation of size scaling and elastic interaction between nano-scale defects in collision cascades.  EPL 110, 36001.  doi:10.1209/0295-5075/110/36001.
  * Yu et al. (2017) Yu, H., Yao, Z., Idrees, Y., Zhang, H.K., Kirk, M.A., Daymond, M.R., 2017.  Accumulation of dislocation loops in the $\alpha$ phase of Zr Excel alloy under heavy ion irradiation.  Journal of Nuclear Materials 491, 232–241.  doi:10.1016/j.jnucmat.2017.04.038.
  * Zarestky (1979) Zarestky, J.L., 1979.  Lattice dynamics of hcp and bcc zirconium.  Ph.D. thesis. Iowa State University.  doi:10.31274/rtd-180813-3520.
  * Zhou et al. (2006) Zhou, Z., Jenkins, M.L., Dudarev, S.L., Sutton, A.P., Kirk, M.A., 2006.  Simulations of weak-beam diffraction contrast images of dislocation loops by the many-beam Howie–Basinski equations.  Philosophical Magazine 86, 4851–4881.  doi:10.1080/14786430600615041.
  * Ziegler et al. (2010) Ziegler, J.F., Ziegler, M.D., Biersack, J.P., 2010.  SRIM - The stopping and range of ions in matter (2010).  Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 268, 1818–1823.  doi:10.1016/J.NIMB.2010.02.091.
  * Zinkle and Was (2013) Zinkle, S., Was, G., 2013.  Materials challenges in nuclear energy.  Acta Materialia 61, 735–758.  doi:10.1016/j.actamat.2012.11.004.