Microstructural complexity and dimensional changes in heavily irradiated zirconium

Using atomistic simulations based on the creation-relaxation algorithm, we explore the evolution of microstructure in irradiated zirconium over a broad range of radiation exposure. In agreement with experimental observations, we ﬁnd that at relatively low temperatures, microstructure evolves towards an asymptotic dynamic steady state forming at doses close to 1 dpa. Simulations show the spontaneous formation of a -type interstitial dislocation loops, gradually transforming into an a -type extended dislocation network and giving rise to macroscopic anisotropic dimensional changes in a textured material. A fully developed a -type interstitial / vacancy and c -type vacancy dislocation microstructure corresponds to the highest degree of irradiation growth anisotropy and vanishingly small volumetric swelling of the material.


I. INTRODUCTION
Zirconium alloys are extensively used in nuclear fission power generation as cladding for uranium oxide fuel. During service, zirconium alloys are exposed to large irradiation doses exceeding 15 displacements per atom (dpa) [1,2]. The dynamic structural transformations stimulated by the highdose irradiation alter the properties of zirconium alloys both at microscopic and macroscopic scales, giving rise to dimensional changes and deformations observed on the scale of entire components.
At operating temperature and pressure conditions, pure α-Zr and its alloys used for cladding adopt a hexagonal close packed (hcp) crystal structure. Let a and c denote the crystallographic directions parallel to 2110 and 0001 , respectively [3]. Claddings are manufactured as highly textured pipes such that a large fraction of grains are oriented with their crystallographic c axes pointing in the radial direction [4], except for CANDU reactors where the c axes point along the circumference of the tube [5]. When subjected to the flux of neutrons generated by the fissile nuclear fuel and with no applied stress, these textured pipes exhibit what is known as "irradiation-induced growth" (IIG): Zr grains undergo a expansion and c contraction, which in the pipe geometry of cladding are referred to as the axial and radial directions, respectively. The experimental data points [6] shown in Fig. 1  that IIG saturates [3]. In the limit of very high dose, typically >5 dpa, and higher temperatures, zirconium exhibits "breakaway growth," where the irradiation-induced strain varies linearly as a function of irradiation exposure. In this work, we focus on exploring the range of radiation exposure and temperatures corresponding to the pre-breakaway mode of microstructural evolution and radiation-induced growth.
There is still no consistent minimum-parameter model for the microscopic picture of evolution of zirconium alloys explaining the fundamentals of IIG across the experimentally relevant range of temperatures and radiation exposure. One of the established and widely used models for IIG was developed by Woo et al. [5,[7][8][9][10], who argued that the observed anisotropy of growth stems from the anisotropy of diffusion of radiation defects. Indeed, in an hcp crystal the diffusion coefficient for thermally activated motion of radiation defects is not a scalar quantity but a nonisotropic tensor. By reformulating the rate theory model for microstructural evolution of isotropic irradiated materials [11], to take into account the anisotropy of thermal diffusion of defects [8], the diffusion anisotropy difference (DAD) model [5,7,8,10] enables making anisotropic predictions for the rates of accumulation of defects in dislocation loops with a-or c-type Burgers vectors.
Still, it is known that the rate theory approximation is valid only in the limit where the concentration of defects is low and temperature is relatively high, imposing a constraint on the range of temperatures where the treatment applies. Using rate theory models, it is difficult to explore microstructural transients and in particular the low-temperature mode of evolution, where microstructure exhibits saturation as a function of irradiation dose [12][13][14]. Furthermore, rate theory models require using multiple input parameters [15], which are often difficult to determine accurately from atomistic simulations [15][16][17][18][19]. At the same time, direct simulations of recoil events initiated by neutrons, and the resulting atomic collision cascades, struggle to reach doses above 1 dpa [20].
FIG. 1. Irradiation-induced growth strain in pure α-Zr in a and c directions as a function of dose. Our simulations predict that the growth strain ε Sim saturates as a function of dose, in qualitative agreement with the variation of experimental strain ε Expt [6] observed at a relatively low irradiation temperature T = 353 K. Simulated growth strain is scaled by a factor of 0.15 for a quantitative comparison. See Sec. II for details on the simulation method.
In this study we use the recently developed creation relaxation algorithm (CRA) to simulate the formation of complex high-irradiation-dose microstructures of hcp zirconium. The algorithm has previously been successfully applied to the investigation of high-dose microstructures of body-centeredcubic (bcc) Fe and W [21][22][23], where in the latter case, a model for a complex stress state of tungsten showed good quantitative agreement with experiment. Below, we apply the method to investigate the formation of complex defect and dislocation configurations in heavily irradiated Zr in the limit of high dose 1 dpa. Figure 1 illustrates how parameter-free CRA simulations predict dimensional changes that exhibit a pattern of variation as a function of dose similar to what is observed in experiment, and shed light on the underlying fundamental processes responsible for the saturation of strain. This mode of microstructural evolution, similarly to iron and tungsten [21][22][23], does not involve diffusion of radiation defects. The evolution of microstructure and dimensional changes occurring at lower irradiation temperatures is driven by the competition between accumulation and relaxation of strongly fluctuating stresses produced by the defects [21][22][23], resulting in a saturated dynamic steady state and asymptotically stationary dimensional changes of the irradiated material [13].
Recently, Tian et al. [24] published a study on irradiated Zr, also employing the CRA, where some of the numerical results are similar to those presented here. There are key differences between our analyses. While Tian et al. provide a useful comparison of the CRA with cascade simulations and explore the details of structural relaxations, we have used the CRA simulations to link the microscopic properties of individual defects to the deformation that they produce in the cladding. Also, we provide an explicit comparison of simulations with experimental observations of deformation and explain why the CRA simulations can be used for exploring phenomena occurring at temperatures pertinent to the operating conditions of a reactor, and not just under the cryogenic conditions considered by Tian et al. Furthermore, we are able to identify an energy barrier for the onset of the breakaway growth. This activation energy is an important parameter defining the temperature and timescales of microstructural evolution of irradiated zirconium.
The paper is organized as follows. In Sec. II, we describe the creation relaxation algorithm and apply it to simulating the evolution of irradiation-induced microstructures at low homologous temperatures. In Secs. III and IV, we analyze the structural properties and complexity of microstructures across a broad range of irradiation doses before summarizing key results in Sec. V.

II. HIGH-DOSE MICROSTRUCTURE
In this section we outline the CRA simulation method, define its range of validity, and compare it with other approaches to modeling microstructures in the high-dose limit. The evolution of microstructure under irradiation in the CRA approach is simulated using a repeating sequence of three steps: (1) Create a small number of Frenkel pair defects by randomly choosing one or more atoms in the simulation cell and placing them to randomly chosen positions elsewhere in the cell. The structure of atomic configurations involved in the execution of the algorithm is not assumed to be crystalline.
(2) Relax the resulting atomic configuration to a local energy minimum.
(3) Repeat the above two steps up to a specified level of damage, each time increasing the cumulative dose φ by the amount N FP /N cell , where N FP is the number of atoms randomly chosen and displaced at a given step of the algorithm, and N cell is the total number of atoms in the simulation cell.
In literature, there are various definitions for the dose φ, which is a measure of exposure of a material to irradiation [16] rather than a parameter characterizing its physical state. A commonly used measure of radiation exposure, or dose, was proposed by Kinchin et al. and Norgett et al. in Refs. [25,26], and this quantity is often referred to as "NRT dose" and expressed in units of NRT dpa, or simply dpa. In the CRA, the number of displaced atoms is known exactly over any number of steps of the algorithm, thus allowing for a natural measure of a "canonical dose," given in units of "canonical displacements per atom" (cdpa) [21]. The canonical and NRT measures of dose can be related by comparing concentrations of defects predicted by the CRA to cascade simulation data [20]. The profiles were found to be almost identical and related by a transformation f (x) → a f (bx) for constant a and b, implying that predictions of physical quantities as a function of canonical dose are transferable to predictions measured against NRT dose through a uniform linear scaling.
Can the evolution of microstructure in the limit of high dose be explored by a direct simulation of recoil events as opposed to using the CRA? To compare the computational efforts involved, consider an example of a collision cascade simulated by molecular dynamics [16]. The average number of defects, or atomic displacements [25], produced in an ideal single crystal by a collision cascade event with damage energy where E d is the threshold displacement energy [26]. E d varies in the interval from 20 to 100 eV and is the average, over orientations and statistical realizations of collision cascade events, energy required to produce a Frenkel pair [27,28]. Since any simulation of a collision cascade involves the eventual dispersion of the energy of atomic recoils into electronic excitations and thermal atomic vibrations, with the characteristic energy of E th ∼ 0.01 eV [29], the minimum number of atoms N cell that a simulation cell has to contain to enable thermalization is N cell T d /E th . Hence, the number of displacements per atom (dpa) produced in a collision cascade event with the damage energy T d is which, according to the above equation, is independent of T d . Assuming E d ∼ 50 eV, we find that a direct molecular dynamics simulation of an event initiated by an atomic recoil produces φ ∼ 10 −4 or less dpa per cascade simulation. Hence, accumulating the dose of 10 dpa directly by molecular dynamics simulations requires modeling in excess of ∼10 5 sequential cascade events. This estimate does not depend on the detailed methodology of molecular dynamics simulations and is independent of the energy of a collision cascade event. This is confirmed by the recent analysis showing that reaching the dose of ∼0.1 dpa required generating ∼1000 sequential cascade events [17,18,20,30]. Performing 10 5 sequential simulations is computationally demanding [17,20], highlighting the fact that simulating high-dose microstructures necessarily requires using the most numerically efficient algorithms possible and making the CRA a method of choice for the initial analysis.

A. Low-vs high-temperature microstructural evolution in α-Zr
Since the total energy of an atomic configuration is minimized at each step of the CRA algorithm, the evolution of microstructure as a function of dose is driven purely by the relaxation of atomic positions, or in other words, by the internal stress resulting from the formation and accumulation of defects. Since in CRA, structural relaxation involves no consideration of the thermal motion of atoms, the approach might appear to involve an element of oversimplification.
However, recent studies performed using the CRA have captured important trends observed in experimental measurements of heavily irradiated materials [21][22][23]. The fundamental reasons for this agreement with observations stem from the difference between the timescales of the relatively slow thermal diffusive relaxation and fast stress-driven atomic dynamics compared to the dose rate; it is in the regime where stress relaxation and not diffusion dominates the microstructural dynamics that CRA provides a good description of accumulation and evolution of radiation damage.
Under the conditions where the dose rate is high for the thermally activated diffusive processes to have a significant effect on the microstructure before new defects are created, yet it is low enough to enable the internal stresses to rapidly relax a newly created atomic configuration into a local energy minimum, CRA simulations are able to provide a reasonably accurate description of the evolving microstructure over a broad range of irradiation exposure, particularly in the limit of very high dose [21][22][23].
Thermally activated spontaneous microstructural transformations occur at a rate τ −1 given by the Arrhenius law [31,32] for the attempt frequency ν, approximately equal to the Debye frequency ν D = ω D /2π , the activation energy for migration E a , temperature T , and the Boltzmann constant k B . Self-interstitial atom defects, hereafter referred to as selfinterstitials, are in general significantly more mobile than vacancies [33,34], and hcp α-Zr is no exception. The activation energy for migration of self-interstitials is as low as 0.008 eV in contrast to a barrier for migration of ∼0.5 eV for vacancy hopping in the basal plane [35]. Even at the operating temperatures of pressurized water reactors (PWRs) in the interval 288 • C to 350 • C [5], Eq. (2) implies that the mobility of self-interstitials is on the order of a thousand times greater than vacancies. The potential barriers for migration of self-interstitials are sufficiently small to be overcome by internal stresses in a material with high-defect content, and this drives the evolution of microstructure in the absence of thermal activation. Vacancies are less sensitive to stress [36] and the potential barriers for migration of vacancies are higher and, thus, in a typical CRA simulation, vacancies do not move appreciable distances unlike the interstitials that migrate and easily cluster. The migration energy of vacancies in hcp zirconium computed using ab initio methods is 0.55 eV in the basal plane and 0.65 eV out of the basal plane [37]. For comparison, the migration energy of vacancies in bcc iron is in a very similar range from 0.55 to 0.65 eV [34,38]. The similarity of vacancy migration energies in Zr and Fe is confirmed by the observed temperatures of stage III of resistivity recovery, which in Zr is close to 280 K [39] and in Fe is close to 250 K [38]. Vacancy formation energies of 2.1 eV in bcc iron [40] and 2.07 eV in hcp zirconium [37] are also very similar.
Most notably, the temperature intervals of thermal recovery of heavily irradiated zirconium alloys and steels also show remarkable similarity. The fully developed high-dose defect microstructure of heavily irradiated zirconium alloys recovers in the temperature interval from 300 • C to 450 • C [41]. Similarly, the mechanical properties of heavily irradiated ferritic steels recover their extreme brittleness over the same interval of temperatures from 300 • C to 450 • C [42,43]. It is this range of temperatures that is of prime interest in relation to CRA simulations since it enables estimating the activation energy for recovery of a high-dose microstructure.
Inverting Eq. (2) and interpreting τ −1 as the thermal recovery timescale, we find the effective activation energy for the thermal self-evolution of microstructure as Using this equation, for T * = 350 • C, τ = 1 h and ω D = 3.67 × 10 13 s −1 [44], we find the activation energy of E a ≈ 2.0 eV. This value is much higher than the migration energy of vacancies in pure crystalline zirconium or pure iron, showing that the recovery of a high-dose microstructure is controlled by higher-energy processes than simple vacancy migration. While highlighting the potentially significant part played by the impurities and solutes, the above estimate suggests that a high-dose microstructure represents a far more thermally stable configuration of an irradiated material than a crystal containing a very low concentration of defects, typically used in resistivity recovery studies [38,39]. The high thermal stability of microstructure in a heavily irradiated material likely stems from the occurrence of vacancy clusters [45].
Using the above estimate for E a , we can define a range of temperatures and dose rates where CRA simulations are expected to apply (see, e.g., Fig. 1 of Ref. [21]). Assuming the dose rate of 10 dpa/year characteristic of reactor operating conditions [13] and the effective activation energy for thermal microstructural relaxation of 2 eV, we find that the stress-driven modes of relaxation dominate microstructural evolution at temperatures below ∼520 K. This condition is well correlated with experimental observations [13,14], showing that the modes of evolution of microstructure below and above this temperature are qualitatively different, and it broadly defines the range of validity of the simulations included in this study.

B. Elastic parameters of defects in α-Zr
The elastic parameters of point defects in Zr can be deduced from the data available in the literature and this can help foresee some simple, but important, features of the CRA analysis presented below. Consider the relaxation volume of a defect , which is equal to the volume difference between a crystal with and without a defect [46]. Relaxation volumes of point defects in metals exhibit a systematic trend, namely, that self-interstitial atom defects cause lattice swelling and vacancies give rise to lattice contraction [47,48]; however, the magnitude of relaxation volume is always larger for the self-interstitials. Thus, the difference in relaxation volumes of self-interstitial and vacancy defects always results in swelling.
The general elastic parameters of a point defect are encoded in the relaxation volume tensor i j , which is a universal quantity characterizing the elastic field of a defect [46,49]. For a defect in a sufficiently large simulation cell of volume V cell , one may compute i j from the strain of the simulation cell ε i j , sometimes called the global strain, associated with full dimensional relaxation of the cell, namely [50] The relaxation volume of the defect equals the trace of this tensor, = ii . Often, authors provide not the relaxation volume tensor but instead calculate the elastic dipole tensor P i j , which is related to i j in the same way as elastic strain is related to elastic stress, namely, where S i jkl is the elastic compliance tensor. Table I contains i j values computed from ab initio data using the elastic TABLE I. Components of relaxation volume tensors and total relaxation volumes computed for vacancies (V) and low formation energy self-interstitials occupying the octahedral (O) and basal octahedral (BO) sites. Due to the symmetry of these point defects, the off-diagonal components of the relaxation volume tensor vanish and 11 = 22 . Data are presented both for the MA#3 potential and ab initio computed elastic dipole tensors [51], transformed using the values of elastic constants from Ref. [52]. Values of i j are given in the units of volume 0 of a Zr atom in a pristine hcp crystal. constants C i jkl from Ref. [52] and P i j for a vacancy and two low-formation-energy self-interstitials [51], together with i j derived from one of the interatomic potentials used in our CRA simulations. The components of S i jkl in terms of C i jkl with hcp symmetry are given in Appendix B.
The character of i j for the point defects is apparent in the relaxed atomic configurations rendered in Fig. 2 where self-interstitials cause significantly larger displacements than vacancies. Interestingly, due to the anisotropic hcp crystal structure, the 33 component of the basal octahedral selfinterstitial is negative, indicating a small negative-c strain. This is reflected in Fig. 2(b) where atoms in noncoincident (0001) planes experience compressive hydrostatic pressure and are displaced towards the self-interstitial by a small amount. Nevertheless, the overall effect is that in both sets of relaxation volumes, | ii | is larger for the self-interstitials than for the vacancies, suggesting that a net swelling is going to occur in CRA as confirmed by the data shown in Fig. 1.
However, in most cases the measurements of IIG indicate near-zero volume change. Experimental data show that in the limit where IIG is particularly well pronounced, vacancies move and cluster, forming dislocation loops. For a sufficiently large dislocation loop, where dislocation core effects may be neglected, the relaxation volume tensor of the loop Dis i j has the form [49] Dis where b i is the Burgers vector of the dislocation forming the perimeter of the loop, and A i is the loop area vector. The corresponding relaxation volume is Dis where we follow a convention for the direction of the Burgers vector b that for a self-interstitial (vacancy) type dislocation loop Dis ii is positive (negative). One may construct a dislocation loop from a collapsed platelet of self-interstitials or vacancies [54].
In the limit of a large loop, for n point defects contained in the loop, the formation volume is equal to ±n (in units of the volume of a Zr atom in the pristine hcp phase) for interstitial (+) and vacancy (−) type loops. For example, this behavior of the formation volume per defect tending to ±1 for large dislocation loops has been confirmed by direct simulations of dislocation loops in bcc tungsten [54,55]. Thus, in a scenario where equal numbers of self-interstitials and vacancies form large dislocation loops, the net volume change is expected to be close to zero.

III. RESULTS AND COMPARISON WITH EXPERIMENT
For all of our CRA simulations, the starting configuration, corresponding to φ = 0 cdpa, is a perfect crystal of hcp α-Zr. Periodic boundary conditions are applied to a hexagonal supercell containing 2 005 056 atoms with edges parallel to We find that the resulting microstructures are similar to those obtained from inserting one Frenkel pair at a time. After each insertion step, the resulting atomic configuration is relaxed to a local energy minimum using the conjugate gradient method. Further details of the simulation cell geometry and convergence testing are given in Appendix A.
Throughout this paper, we employ the three embedded atom method potentials developed in Ref. [53]. In correspondence to the numbering provided by the authors of the potentials, we refer to each potential as MA#1, MA#2, and MA#3. The potentials were obtained from the NIST Interatomic Potentials Repository website [56][57][58], noting that we have employed the 24 Feb 2009 updates to MA#2 and MA#3 that were designed to better treat the radial distances smaller than 0.5 Å [56]. A survey of the literature indicates the MA#3 potential appears to be the most frequently used for simulating irradiation damage [18,19,30,[59][60][61][62][63][64][65]. Commonly cited reasons are that MA#3 was fitted to ab initio point defect formation and stacking fault energies. Nonetheless, some authors make use of MA#2 [66]. Unless otherwise stated, all the figures shown in the main text are produced by the MA#3 potential, with data for the MA#1 and MA#2 potentials presented in the Appendixes.
The LAMMPS [67] molecular dynamics program was used to construct and relax the atomic configurations [68]. Data were visualized and processed with the OVITO [69] and PAR-AVIEW [70] software.

A. Saturation of global stress and strain
The shape and volume of the supercell remains fixed throughout the simulation. As such, as the dose increases, the supercell acquires a global internal stressσ i j . The diagonal components of internal stress, plotted over 10 cdpa in Fig. 3(a), are always expansive and increase rapidly from 0 to 2 cdpa before saturating, while the off-diagonal components are negligibly small in comparison. The expansive character of internal stress is consistent with the fact that an application of external pressure is required to keep the boundaries of the cell fixed, and the components of this external pressure are all positive [71]. Interestingly, while the insertion of Frenkel pairs must break the symmetry of the pristine configuration, over all canonical dosesσ i j approximately retains the form of a rank-2 material tensor with axial symmetry such thatσ 33 is distinct fromσ 11 ≈σ 22 . The prismatic plane normal stresses are initially approximately equal and diverge slightly at 3 cdpa before attaining similar values at canonical doses >5 cdpa. On the other hand, theσ 33 is significantly smaller and exhibits relatively smaller statistical fluctuations.
If we were to allow the cell shape to relax freely, in equivalence to using traction-free boundary conditions, then the supercell would acquire a global strainε i j . As the global stresses are small relative to the elastic constants of α-Zr, we may findε i j using Hooke's law such that whereσ kl is the global internal stress developing in a cell with fixed boundaries. The validity of using Eq. (7) to compute cell strain instead of relaxing the cell under zero stress conditions has been confirmed in Ref. [24]. The resulting strains are plotted in Fig. 3(b) where, as observed in experiment [6], saturation occurs at high dose and a textured material exhibits expansion in the basal plane and contraction in the c direction. As discussed in Appendix C, for small canonical dose <0.1 cdpa, the exact shape of the strain-dose profile curve is potential dependent.
The strain versus dose profiles found in our simulations are similar to the measurements in the pre-breakaway growth regime of IIG in single crystal α-Zr [6,13]. In particular, the character in strain exhibiting a expansion and c contraction is consistent, and agrees with observations. Furthermore, the apparent saturation in strain is also observed at doses >1 dpa before the breakaway growth occurs. Electron microscope observations of radiation defects produced in zirconium by ion irradiation also exhibit saturation at doses above 0.3 dpa at 300 • C [14]. In the next section, we describe the underlying microstructure and relate it to the global stress and strain parameters.

B. Steady-state microstructure
In a CRA simulation, atoms are rearranged via the gradual accumulation of Frenkel pair defects. These point defects may annihilate the point defects already present in the microstructure, stabilize in isolation, or cluster into larger defects. A notable type of defects arising from clustering of selfinterstitial point defects are interstitial-type dislocation loops. As discussed in Sec. II, the predominant absence of vacancytype dislocation loops and clusters is a consequence of the low-vacancy mobility pertinent to the CRA.
Point defects are identified using the Wigner-Seitz method. For all the three potentials employed in simulations, vacancies are distributed approximately uniformly in the simulation cell across all the canonical doses. On the other hand, a majority of the highly mobile interstitials cluster, to form dislocation loops or full atomic planes. When visualizing the defect structures, we therefore remove atoms that are both reported as interstitials by the Wigner-Seitz analysis and identified by an Ackland-Jones structure analysis [72] to be in an hcp environment. This analysis yields the profile in Fig. 4(a) where, in the dilute limit for very low canonical doses <0.1 cdpa, the self-interstitial defects and vacancies stabilize in isolation. While vacancies still tend not to cluster at higher doses, isolated self-interstitials reach a maximum concentration at ∼0.1 cdpa at which point they start clustering into dislocation loops, attaining a stable concentration at >1 cdpa.
Dislocations were identified by the DXA algorithm [73]. The dislocation density evolves through three distinct regimes as a function of the canonical dose, as illustrated by the renders in Fig. 5. For 0.1 cdpa, interstitials coalesce into small but numerous dislocation loops. In the domain of >1 cdpa, the small loops coalesce into larger dislocations forming a system-spanning dislocation network. At even higher doses >2.5 cdpa, this dislocation network equilibrates with the arriving self-interstitials and partially breaks up into a dynamic steady-state microstructure involving isolated defects, dislocation loops, and extended dislocations. Microstructure renders for the MA#1 and MA#2 potentials are shown in Figs. 6 and 7, respectively. Across all the potentials and doses we find that a significant proportion of dislocations (∼70%) are of a interstitial type with Burgers vectors of 1 3 1210 (perfect) and 1 3 1100 (partial). Within the pre-breakaway growth regime at <5 dpa [5], a interstitial-type loops are also observed in experiment [10]. The remaining proportion of dislocations (∼20%) detected by the DXA algorithm have unusual Burgers vectors, e.g., 1 18 4 043 . The dislocations either form closed loops or networks or terminate at another defect, often at twin boundaries. Occasionally, the DXA identifies dense networks of small dislocations on the twin boundary surfaces which leads to sharp fluctuations of the dislocation density profile in Fig. 4(b). The volume fraction of these twinned regions is plotted in Fig. 4(c) for the MA#3 potential where the profile consists of a peak at low canonical dose before saturating at larger doses. The atomic structure of the most frequently occurring type of twin is shown in Fig. 8, which we have identified as a {1121} tensile twin generated by a misorientation of ∼34.84 • about 1010 . Figure 9 contains plots of the dislocation density and volume fraction of twinned regions for all three potentials, where it is clear that twinning is particularly prevalent in the MA#2 and MA#3 potentials but rare in the MA#1 potential. Moreover, the peak in MA#3 twinning is absent in the other potentials which instead display a peak in dislocation density, much alike the profile calculated with the CRA for bcc Fe [21]. For the MA#1 and MA#2 potentials, small a-type dislocations form in large numbers before coalescing into larger loops and networks, giving rise to the peak in dislocation density. However, for the MA#3 potential small dislocation loops with unusual Burgers vectors instead coalesce into twinned regions, suppressing the initial rise in dislocation density and giving rise to the peak in twinning. Although we note that the same types of twins are observed experimentally in Zr alloys [74], twinning is not reported in experimental studies of irradiated zirconium alloys, and thus we believe that their presence might be an artifact of the potentials.
As plotted in Fig. 4, interstitial and vacancy concentrations, dislocation line densities, and the volume fraction of twinned regions saturate at dose φ 1 cdpa. While the defects are continually produced at each step of a CRA simulation, this flux of defects is eventually balanced by the stress-driven recombination together with an increasing probability that at least one of the interstitial and vacancy in a given created Frenkel pair undergoes immediate recombination. Finally, it is important to appreciate the strongly spatially fluctuating distribution of stress underlying this highly damaged microstructure in the material. While the average global stressσ i j is shown in Fig. 3 and is a gradually varying function of dose φ, Fig. 10 illustrates the strongly varying nature of stress σ i j (x) as a function of position x in a cross section of the 5.00 cdpa MA#3 microstructure (cf. the bottom of Fig. 5). Here, we have drawn a heat map (left) of the von Mises stress σ vM (x) which is an invariant quantity given by [40] where repeated indices imply summation. The surface map (right) represents σ vM (x) as a function of position within a two-dimensional cross section of microstructure of a heavily irradiated zirconium. The sharp peaks of large σ vM (x) make the strongly fluctuating nature of the stress manifest. This landscape of stress implies a strongly fluctuating distribution of strain in the material, which is critical to the interpretation of x-ray diffraction measurements of heavily irradiated Zr and Zr alloys.

IV. DISCUSSION
As noted in Sec. I, the trend exhibited by the simulated strain vs canonical dose profile is similar to growth measured  in experiment at relatively low irradiation temperatures, where the growth strains saturate at high dose >1 dpa and exhibit positive-a and negative-c strains. This behavior is correlated with the saturation of defect populations, a key component of which is the production of interstitial-type dislocations with a-type Burgers vectors, as is also observed in experiment. This broadly falls in line with our discussion in Sec. II A, where we argued that the CRA does provide a suitable description of irradiation damage at temperatures below ∼520 K with dose rates of ∼10 dpa/year.
Although the shape of the growth curves is consistent, there is a difference in scale between experiment and the strains predicted by the CRA simulations. This difference was also observed in other studies [22], where it was found that apart from the overall scaling, the observed variation of the predicted strain versus dose was in agreement with observations. Furthermore, the small-c strain relative to a strain gives rise to swelling unlike the near-zero change in volume measured in IIG. The positive swelling predicted by CRA simulations is due to the lack of vacancy mobility and the asymmetry in dimensional changes produced by interstitials and vacancies. This is in line with expectations laid out in Sec. II. Curiously, we also note that for some interstitials, the 33 component is slightly negative (e.g., for the basal octahedral interstitial in Table I). This is true for both ab initio data and potentials used in our simulations. Thus, this presents another mechanism through which the small-c strain appears in our simulation, the finer details of which will be important to explore in future work.
The difference in magnitude between the CRA and experiment strains arises primarily from the partial thermal relaxation of microstructure, present in a real material [22,23]. Indeed, upon annealing a 1.57-cdpa microstructure at 250 • C with fixed cell shape for 1 ns, we found a 67% decrease in the dilatationε ii . The absence of breakaway growth in the CRA simulations is likely due to the lack of temperature effects, combined with the persistent generation of defects. It is well known that breakaway growth is strongly correlated with the emergence of vacancy-type c loops [10]. Thus, given the low effective vacancy mobility and therefore low susceptibility of vacancies to coalesce into dislocations loops in our CRA simulations, the emergence of a significant number of vacancy loops will likely feature only over a large interval of dose, not yet accessible to a direct CRA analysis.
In order to shed light on the role of thermally activated processes, we estimate the effect of vacancy-loop coalescence on the strains plotted in Fig. 3 as follows. Consider the effect of rearranging vacancy point defects into large dislocation loops. With reference to Table I, for the MA#3 potential the relaxation volume tensor of a vacancy is very close to an isotropic form i j = 1 3 V δ i j , with the relaxation volume V = −0.38 0 , where 0 is volume of an atom in an hcp lattice. Making use of Eq. (6) and the experimental observation that c loops (b i ∝ δ i3 ) are observed to be vacancy type with the habit plane parallel to the basal plane (A j ∝ δ j3 ), we may write the relaxation volume tensor of a c-vacancy loop as i j = −nδ i3 δ j3 , where n is the number of vacancies in the loop. The effect of a rearrangement of individual vacancies into vacancy-type c loops on global strainε i j may be estimated using the equation for vacancy concentration c V , given in the units of vacancy per lattice site, and the fraction of vacancies f that coalesce into c-type loops. By using the vacancy profile shown in Fig. 4(a) for c V , with f = 0.1 one may recover a profile of strain as a function of dose that is similar to that observed in experiment as shown in Fig. 11. Thus, we prove that it is the lack of vacancy mobility in the low-temperature conditions that is responsible for the relatively small c strains and irradiationinduced dimensional changes predicted by simulations.

V. CONCLUSIONS
In summary, in this study we have explored how the microstructure of irradiated hcp α-Zr evolves over a broad range of exposure of the material to irradiation. We employed the CRA method that both affords computational efficiency and requires minimal parametrization to simulate microstructures corresponding to high dose. In the CRA, the relaxation of atomic positions is driven exclusively by internal stresses and, thus, it is striking that with no diffusion of defects included in the treatment, our simulations yield global and microstructural properties that are similar to those observed in experiment, including the following: (i) The natural occurrence of positive-a and negative-c strains.
(ii) Dense defect microstructures containing dislocations with a-type Burgers vectors.
(iii) Saturation in strain and defect populations at high dose.
The microstructure at all doses is found to be characterized by a high degree of complexity, and its evolution is driven by the strongly fluctuating internal stresses.
Our work also offers clarity about the origin of the transition between the observed low-and high-temperature modes of microstructural evolution [13,14], distinguished by the relative timescales of dose rate, and stress-driven and thermally activated processes. This classification enables identifying FIG. 11. An estimation of the effect of vacancy-type loops on c strain. 10% of vacancies were assumed to rearrange into large-c vacancy-type loops, thus altering the strains presented in Fig. 3 according to Eq. (9). further points that require understanding for the full treatment of underlying mechanisms responsible for the breakaway growth and the role of alloying elements in delaying such rapid increase in radiation-induced internal strain. While it is clear that the migration of vacancies represents an essential element of evolution, stimulating the breakaway growth and vacancy loop formation, specifically how this mechanism leads to the observed behavior at high temperature remains an open question. As we discuss in Sec. II, the ∼2-eV energy implied by the temperature interval separating the lowand high-temperature regimes, is significantly higher than the ∼0.5-eV migration energy of an isolated vacancy. The question about the origin of this high apparent activation energy remains outstanding. Experimental observations strongly suggest this activation energy is a fundamental characteristic parameter of a heavily damaged material and hence the question warrants extensive further investigation.
To obtain further information on the data and models underlying the materials described in the paper, please see Ref. [75]. and ( 2 3 , 1 3 , 3 4 ). When repeating the primitive cell, the number of repetitions along a 1 is equal to that along a 2 . Let N a and N c denote the number of repetitions along a 1,2 and c, respectively. Due to the periodic boundary conditions, in order to ensure that the distance between a given atom and its periodic image in a neighboring supercell is approximately equal along each cell vector, we require that |a|N a = |c|N c (for |a| := |a 1 |). Thus, for a given choice of N a , we determine N c such that where Round(x) rounds x to the nearest integer. In practice, we used |a|/|c| ≈ 3 8 for all the three potentials which, for increasing cell size, results in negligible loss of accuracy and enables a comparison between different potentials using the same number of atoms. For a given choice of repeating the primitive unit cell, we denote the resulting supercell N a × N a × N c .
In Fig. 12, we show how the profile of pressure in the supercell [ Fig. 12(a)] and the change in energy per atom [ Fig. 12(b)], where changes are calculated with respect to the pristine hcp crystal as a function of canonical dose vary with increasing supercell size. The supercells 12 × 12 × 7, 25 × 25 × 16, 118 × 118 × 72, and 201 × 201 × 123 correspond to 2000, 20 000, 2 005 056, and 9 938 646 atoms, respectively. The fluctuations in these quantities are significantly reduced at the 118 × 118 × 72 supercell size employed in this simulation. Furthermore, while as large a supercell as possible is desirable, with the computational resources available to us, this particular size allowed for the CRA simulation to reach high doses within a reasonable amount of time.
Next, we consider the effect of inserting 10 −3 cdpa per creation step which, for the 118 × 118 × 72 supercell, corresponds to 2005 Frenkel pair insertions. Figure 13 shows strains as a function of canonical doses for one Frenkel pair and 10 −3 cdpa defect insertions per creation step of the simulation. We note that the profiles are similar for both choices. Furthermore, we find that the resulting microstructures are also similar, as was found in recent studies of irradiated bcc tungsten [22,23]. Thus, in favor of computational efficiency, we insert multiple Frenkel pairs per creation step.
The elastic compliance tensor S i jkl is related to C i jkl via S i j pq C pqkl = 1 2 (δ il δ jk + δ ik δ jl ), and, due to also being a fourth-rank tensor with the same symmetries, must have the same tensorial form as that given for C i jkl by Eq. (B1). This may be written explicitly by substituting the C i jkl components in Eq. (B2) with the corresponding S i jkl such that where ζ = C 3333 (C 1111 + C 1122 ) − 2C 2 1133 .

APPENDIX C: SIMULATION DATA FOR DIFFERENT POTENTIALS
Here we present and compare the results of applying CRA to Zr using all the three MA#1, MA#2, and MA#3 potentials. It is apparent from Figs. 14 and 15 that saturation occurs for all three potentials and despite the differences in magnitude, the stress and strain profiles are similar. MA#3 produces smaller strains than the other two potentials yet has a similar pointdefect concentration profile.
Key differences between the potentials lie in the profile of dislocation density and formation of twinned regions, as discussed in Sec. III, and the behavior of c strain at very low canonical dose. At damage levels 0.03 cdpa and 0.07 cdpa for the MA#2 and MA#3 potentials, respectively, a small positive strain accumulates in the [0001] direction, acquiring a maximum value of ∼0.1% at ∼0.01 cdpa before decreasing and changing sign at ∼0.06 cdpa. This is in contrast to the immediate c contraction found for the MA#1 potential. Figure 16 shows that, even with creating only one Frenkel pair per insertion step, this behavior remains unchanged. In this regime where only isolated point defects populate the simulation cell, the difference between the potential originates from the difference in the point-defect formation energies and deformation properties.