Dislocation density transients and saturation in irradiated zirconium

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$.


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).
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 MeV 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°C 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 Highest dose in each sample and highest dose rate (12 × 10 5 dpa/s) Lowest dose in each sample and lowest dose rate (0.7 × 10 5 dpa/s) 0.1 × CRA 0.15 dpa XRD 0.50 dpa XRD 1.00 dpa XRD 2.00 dpa XRD 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 -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).
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 ∼15 dpa 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 ( -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°C to 350°C, exhibiting an hcp-bcc instability only at significantly higher temperatures above 860°C (Willaime and Massobrio, 1989).
With increasing dose, the interstitial and vacancy-type dislocation loops with the 1 3 ⟨2110⟩ 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' ⟨2110⟩ and contraction along the 'c' ⟨0001⟩ crystallographic directions, saturating at doses larger than 1 dpa (Holt, 1988). At temperatures above ∼300°C and doses exceeding ∼5 dpa, large vacancy-type 'c loops' with Burgers vectors 1 2 ⟨0001⟩ or 1 6 ⟨2023⟩ 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 -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 dpa before saturating at larger doses. The microstructures produced by our simulations indicate that at very low doses ≪1 dpa, small dislocation loops form, which subsequently grow and coalesce into a complex interconnected dislocation network developing at moderate doses <1 dpa. 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 ≫1 dpa. 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 selfinterstitial cluster population size distribution follows a power law ( ) ∝ 1∕ 2.2 developing in the saturation regime corresponding to high neutron exposure. Below, we show that this fact, also confirmed by experimental observations , 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.

Dislocation density
There are multiple ways to define the density of dislocation lines in a deformed or irradiated material. For all the measurements and computations in this study, is defined as a scalar ratio of the total dislocation line length to volume containing the dislocations, namely (Hull and Bacon, 2011) Here is a position vector on a dislocation line such that d = d for unit tangent vector and arc length , and the integration is performed with respect to the arc length over all the dislocation lines ⟂ in . 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 that measures the number of dislocation lines crossing an open surface as (Hull and Bacon, 2011) where d is a vector area element of 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 of the dislocation. Significant physical quantities such as dislocation energy and the Peach-Koehler force both depend on . The Nye tensor is a tensorial measure of dislocation density that is a linear function of lattice curvature (Nye, 1953) and is also a function of position such that (Jones et al., 2016) where ⊗ denotes the tensor product, integration is performed over all of the dislocation lines in the system, and ( ) is the Dirac delta distribution defined by the property for an arbitrary well-behaved function ( ). Whilst full information about the dislocation content is contained in Equation 3, attempting to average over a volume can be problematic. Essentially, this stems from the fact that the integral of d along a dislocation segment contained in that starts and ends at and respectively is − . 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 microbeam 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 in Equation 1 which is the same definition as that employed in our X-ray line profile analysis.

Experiment
Dislocation densities were measured in four 3 mm × 1 mm × 0.5 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 × 10 14 m −2 and a characteristic 'split-basal' texture due to processing, where the basal poles are aligned along the normal direction (ND), with a ±30 degree tilt towards the transverse direction. Irradiation induced growth strains in similarly textured Zircaloys are known to saturate at doses below ∼10 dpa at 320°C (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 MeV protons at 350°C 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 ±10°C 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 eV and 40 eV respectively (Stoller et al., 2013). A typical dose vs. depth profile in one of our Zircaloy-4 samples consists of a plateau region extending ∼10 µ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 ∼30 µ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 dpa respectively. Within the first 30 µ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 µm × 100 µm cross section, the samples were scanned in cross-section from the surface to a depth of 50 µm within the sample at 2 µ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 µ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.
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 ( ), where 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 LaB 6 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 , denoted , are related to the strain distribution by where is the Fourier variable and ⟨ 2 , ⟩ is the mean-square strain. Wilkens (1970) derived an expression for ⟨ 2 , ⟩ 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 perpendicular to the line direction. A number ⟂ 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 . The dislocation density 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 = √ . Whilst is one of the fitting parameters in the CMWP software, thus affecting the value of , this article is principally concerned with the determination of 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) where is a parameter that is refined by the profile fitting CMWP algorithm, termed the dislocation contrast factor. The value of was evaluated for dislocations in -Zr by Balogh et al. (2016). The Wilkens function ( ), where = ∕ , has the following asymptotic forms in the limit of small and large , respectively The numerically obtained formula for ( ) 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 Topping et al., 2018) as well as other microstructural features  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 nm −1 < < 13 nm −1 were included in the fitting procedure and the uncertainty in 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 2b 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.

Simulation
The accumulation of defects and microstructural evolution as a function of dose was simulated using the Creation Relaxation Algorithm (CRA) . 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 -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 -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 ∼300°C, see Adamson et al. (2019). At temperatures above ∼300°C, 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 ∼300°C in proton irradiated Zircaloy-2 when samples irradiated to 2 dpa were annealed for 1 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 ∼300°C.
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 for the processes primarily responsible for the observed thermally activated behaviour must be close to ∼2 eV. Also, as described in §2.2, the proton-irradiation defect production dose ratė 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̃ 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 , using the dose rate model by Nordlund et al. (2018), we find that the two rates are comparable iḟ where is the threshold displacement energy required for forming a defect, the attempt frequency is ≈ ∕2 = 5.84 × 10 12 s −1 , given the Debye frequency = 3.67 × 10 13 s −1 (Zarestky, 1979) and = 0.861 × 10 −4 eV/K is the Boltzmann constant. Taking = 2 eV and = 40 eV, and solving equation (8) for̃ , we find̃ = 625 K≈350°C. 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̃ , the defect structures generated by irradiation evolve predominantly through fast athermal stress relaxation .
The above esitmate for̃ ∼350°C 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 ∼300°C 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 -Zr (Varvenne et al., 2014) of ∼0.5 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°C.
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  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 Boleininger et al., 2022). The latter forms spontaneously at doses above approximately 0.3 dpa (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 ≈ 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 selfinterstitial 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 . 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  to the high-temperature mode dominated by the Arrhenius thermally activated diffusion (Allnatt and Lidiard, 1993), in zirconium alloys spans approximately from 300°C to 400°C, 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). Selfinterstitial 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°C to 400°C. Over this temperature interval, the mode of microstructural evolution changes from that dominated by stress fluctuations and coalescence of selfinterstitial 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) 1 . 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 ∼ 2M and ∼ 10M atoms, with the cell edges parallel to the [2110], [1210] 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 meVÅ −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 ∼1 GPa, 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 for wavevector magnitude , and the sum runs over all the pairs of atoms positioned at and separated by distance = | − |. 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 for a system of atoms scales unfavourably as 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.
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).

Dislocation structure and distribution
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 (10) vacancies that are not large enough to relax into dislocation loops. We do find a small number of much larger clusters containing (100) Figure 4: Simulated {1102} 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 ∼ 10 15 m −2 .
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  3a), small loops form before coalescing into a dislocation network (Figure 3b) at which point the dislocation density saturates. At high doses (Figure 3c), 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.
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(Mason et al., , 2021Warwick 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. (2020Mason et al. ( , 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 {1102} 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 . 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.
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 5a shows the frequency distribution of clusters containing interstitials as calculated by OVITO for the doses corresponding to the renders in Figure 3. For clusters with population sizes < 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, ) 5 < N < 10 3 100 < N < 10 3 Figure 6: Estimation of dislocation density from interstitial cluster perimeter density.
corresponds to the number of atoms in the cluster and the perimeter of a given cluster is estimated according to Equation 13. simulations employing the CRA have shown evidence for self-organised critical behaviour  and thus it is reasonable to test the cluster population size power law distribution hypothesis such that the probability mass function for clusters containing interstitials is given by where the Riemann zeta function (Heynsworth and Goldberg, 1965) is defined as ( ) = ∑ ∞ =1 1∕ for exponent > 1. We may define an exponent that best fits the data shown in Figure 5a via a maximum likelihood estimation. The likelihood function returns the probability of observing the tot measured data points { | ∈ [1.. tot ]} if they were produced from a distribution ( | ) of given parameter(s) . The MLE is produced by determining the value of MLE maximising ln  such that: Equation 12 was solved numerically using the python package POWERLAW (Alstott et al., 2014) and as shown in Figure  5b, the MLEs and associated standard error for 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 across all three MA potentials. Averaging across the potentials for doses larger then 1 cdpa in the saturation regime, we find the MLE for exponent = 2.28 ± 0.01 for a 2M simulation cell and = 2.23 ± 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). 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 lying in a range min < < max correspond to a circular platelet of interstitial atoms i.e. a dislocation loop (Gilbert et al., 2008). The min and 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 atoms with volume Ω equal to the atomic volume of -Zr. Using the formula ( ⋅ ) for the volume of a dislocation loop, its perimeter can then be estimated as where is the length of the Burgers vector along the normal to the loop habit plane. Here, we assume that the platelets form -type dislocation loops such that is equal to the lattice parameter. Summing the perimeters corresponding to cluster population size distributions between various choices of min and 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 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 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 < < 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 where it is apparent that, assuming that all the loops have the same Burgers vector , ∝ 2 . Treating as a continuous variable, we derive the probability density function as a function of loop diameter where is the exponent entering Equation 10. Hence the diameter of circular loops is expected to be power law distributed with exponent = 2 −1. Using the data derived from the CRA simulations in the high dose limit illustrated in Figure 5b for a 10M atom cell, we find that = 3.46. In experiments performed by  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 ≤ ≤ 4, which agrees well with the values derived from our simulations. Concluding this section, we note that power law exponent values = 2 and = 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 < 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 towards and above the limiting value of 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 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 containing atoms with total energy ( ) we may measure the excess (stored) energy as where ℎ is the cohesive energy of -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 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 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 1 2 0 where Hooke's law determines the stress = 0 and the elastic constants tensor is . Enforcing zero global strain requires that the integral of the strain over the body with volume is zero or, equivalently, that one has zero volume averaged strain This boundary condition is satisfied by the strain Applying Hooke's law to Equation 17, we observe that the system is under a state of global stress ⟨ ⟩ = ⟨ 0 ⟩ and we may treat this as the stress that develops in our simulations. Therefore, we can compute the stored elastic energy where is the elastic compliance tensor (Warwick et al., 2021) related to by = 1 2 + . We find across all the employed potentials that 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 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.
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 where, for dose , the number of vacancies is denoted and 1 is the formation energy of a single vacancy.
Values of 1 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 vacancies identified by Wigner-Seitz analysis and subsequently relaxed an -atom supercell of pristine -Zr containing the same number of vacancies in the same positions. The simulation cell was relaxed and the formation energy was computed as where is the resulting total energy. Figure 7a 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 is due to isolated vacancies. The contribution to from the formation energy of small interstitial clusters and dislocation content may now be calculated as In Figure 7b 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 is shown in Figure 7c 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  Figure 7: Contributions of microstructural defects to the excess energy. Results shown are for the MA1 potential using 10 million atoms. 7a: Comparison of vacancy energy estimated by assuming all vacancies are isolated with formation energy 1 (blue curve) with calculated via explicitly relaxing pristine -Zr containing the same number and arrangement of vacancies. 7b : Breakdown of total excess energy into elastic , interstitial and vacancy contributions. 7c: 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 cdpa.
indicated stored energies associated with irradiation induced defects to be on the order of 0.1 Jg −1 . Their analysis also allowed the authors to infer the presence of defects that are invisible to TEM imaging. At 0.07 cdpa in Figure 7c, we find that the peak in specific energy associated with is 25 Jg −1 that subsequently drops and plateaus at high dose to 12 Jg −1 . At elevated temperatures, the defect content is likely to be ∼ 10% of that calculated in our simulations cf. (Mason et al., 2020(Mason et al., , 2021 and thus we expect the associated stored energy in Zr to be comparable to 1 Jg −1 .

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 dpa. 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 ∼1 dpa, 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 ≈ 2.2. As a function of defect diameter, this results in a power law distribution of defect clusters with exponent of ≈ 3.5, which compares favourably with the range of values 3 ≤ ≤ 4 derived from experimental observations .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 ∼ 350 • C is clear, the dislocation density saturates as a function of dose.

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).