Athermal swelling and creep of heavily irradiated iron under uniaxial stress

The athermal irradiation-induced swelling and creep in iron under the inﬂuence of external uniaxial stress were investigated through atomic scale simulations using the creation-relaxation algorithm. The defect relaxation volume density tensors (or eigenstrains) were evaluated as a function of external uniaxial stress. When the dose exceeded the level at which isolated defects were formed, interstitial-type defect clusters were formed with polarization, resulting in crystal growth in the direction where tensile stress was applied, and growth in the other two perpendicular directions when compressive stress was applied. The concentration of vacancies, isolated self-interstitial atoms, Laves phases clusters, and dislocations in the microstructure were largely unaffected by external stress up to ± 1 GPa. Biased crystal growth was primarily attributed to the anisotropic formation of new lattice planes through the coalescence of interstitial defect clusters, leading to plastic deformation depending on the direction and magnitude of the external stress.


I. INTRODUCTION
Iron and steel are versatile materials with numerous applications in both domestic and industrial settings. Reduced activation ferritic-martensitic (RAFM) steels are chosen as the structural materials of advanced fission and fusion reactors [1,2]. RAFM steels have a body-centred cubic (bcc) crystal structure, which is the same as pure iron. A main reason for choosing RAFM steels is because they show less swelling under neutron irradiation compared to austenitic steels with face-centered cubic (fcc) structure [3]. Commercial ferritic steels are highly resistant to swelling even at extreme exposure to irradiation, showing less than 2% volume change under irradiation at 420 • C up to 200 DPA [4].
Inside a fission or fusion reactor, nuclear reactions generate high-energy neutrons which can penetrate deep inside reactor components. These neutrons are scattered by atoms constituting the component, during which they transfer a fraction of their kinetic energy to the recoiling atoms. If sufficient energy is transferred, a recoil atom becomes a ballistic particle in itself, initiating a cascade of atomic collisions. The total exposure of a material to irradiation is measured in units of displacement per atom (DPA) [5], which is the number density of atoms ballistically displaced by collision cascades. In metals, accumulation of atomic displacements results in the formation of complex microstructures containing defects, such as self-interstitial atoms (SIA), vacancies, voids, dislocation loops, and dislocation lines [6][7][8][9][10][11][12]. Evolution of microstructure can alter both the thermal [13][14][15] and mechanical [16,17] properties of materials, shortening the lifetime of components.
Neutron irradiation experiments of RAFM steels up to 70 DPA in a fission reactor at 330-340 • C [18] show irradiation causing hardening and embrittlement, with saturation of these effects occurring at high doses. Dislocation loops with Burgers vector of a 2 111 and a 100 are commonly observed [4,19], with their relative abundances depending on irradiation condition and alloy composition. RAFM steel will be used as the first wall material of the demonstration fusion plants, known as DEMO, which are expected to experience lifetime neutron irradiation doses between 1 to 10 DPA [20]. Also, Gilbert et al. [21] performed neutron transport calculations. Depending on the design and position of the first wall made of Eurofer, the rate of damage accumulation can range from 8 to 20 DPA per full power year.
In real operating conditions, components of fission and fusion plants can be subjected to applied tensile or compressive stresses due to gravity or magnetic fields. Misalignment of components due to human operation or swelling of materials in limited space can also cause enormous stress to materials, leading to stress-induced anisotropic dimensional changes referred to as creep [22]. Numerous phenomenological models [23,24] based on assumed microstructural evolution were proposed to explain irradiation-induced creep and swelling.
As it is technologically challenging to monitor changes in materials' properties in an extreme radiation environment in situ, computer simulations have the potential to improve our understanding of materials' behavior under irradiation. While state-of-the-art molecular dynamics (MD) simulations of consecutive collision cascades can feasibly reach doses of 1 DPA [25,26], these types of simulations are computationally intensive. For system sizes that can accommodate the formation of extended dislocation microstructure, requiring to the order of a million atoms [27], it is prohibitively expensive to conduct a comprehensive collision cascade study across a parameter space with multiple stochastically independent simulations.
Furthermore, although MD collision cascade simulations, in principle, enable studying thermal-activated relaxation processes, the simulated time (∼ns -μs) is typically at least eight orders of magnitude shorter than experimental time scales (∼min -year). Hence it remains a challenge to study temperature and dose-rate effects using this method. A notable exception occurs in cases where the thermal diffusivity of defects is suppressed due to factors such as low temperature and impurities; in this athermal limit, evolution of irradiated microstructure is driven by the competing effects of defect production and defect recombination, with recombination occurring either through the mutual elastic interaction of defects, or in the case of larger primary recoil energies, by melting and recrystallization of the cascade heat spike region [28].
Recent atomic scale simulations adopted the creationrelaxation algorithm (CRA) [29], which is an accelerated simulation method for simulating irradiated microstructure specifically in the athermal limit. In CRA, an atomic displacement is generated by choosing a random atom and moving it to an arbitrary position, generating a Frenkel pair in the process and incrementing the dose by 1/N DPA, where N is the number of atoms in the system. Subsequently, the atomic configuration is relaxed through minimization of the total energy. CRA presents a simplified picture of the generation of Frenkel pairs by low recoil energy collision events in the athermal limit, achieving a considerate speed-up over conventional collision cascade simulations by replacing explicit time-propagation of recoil atoms with successive displacement and energy minimization steps. We note that simulations using collision cascades, CRA, or a mixture of both methods, have shown qualitative and even quantitative agreement with experiments in the athermal limit [28,[30][31][32][33].
Using CRA, Derlet and Dudarev [29] studied the saturation of defects production in Fe and W. They also studied the spatial fluctuation of stresses and the spontaneous reorganization of microstructure. Tian et al. [34] and Warwick et al. [35] studied the anisotropic swelling of zirconium due to the formation of a-and c-type dislocation loops. Mason et al. [31] explained the negative lattice strain being detected by spatially resolved x-ray Laue diffraction in self-ion irradiated tungsten. Using a combination of CRA and MD, Mason et al. [32] accurately estimated the deuterium retention capacity in heavily irradiated tungsten. Chartier and Marinica [36] used the Frenkel pair accumulation (FPA) method to study highdose irradiation damage in bcc iron. FPA is similar to CRA, but it relaxes the atomic configuration through dynamic evolution. They observed nucleation of C15 clusters at an early stage, which can transform into a 2 111 and a 100 dislocation loops.
Recently, Reali et al. [37] suggested that one can compute the macroscopic stress and strain on the component scale using the defect relaxation volume density tensor ω, which can be obtained from atomic scale simulations, such as CRA or MD. They proved analytically that the relaxation volume density tensor is equivalent to the eigenstrain * of defects [38]. Since the spatially varying eigenstrain can be interpreted as a source of an effective body force [37], one can solve the constitutive equation, and calculate the dimensional changes, through the finite element method (FEM).
In this work, we adopted the CRA to study the effect of external uniaxial stress on defect production and crystal growth in pure iron under athermal irradiation conditions that is realized as a particular form of irradiation-induced swelling and creep. We chose iron as it is the major component of steels and has the same crystal structure as RAFM steels. We will examine the effect quantitatively by calculating the defect relaxation volume density tensor ω (or eigenstrain * ). We will also study the microstructural changes in detail. The anisotropic changes of ω will be explained through observation of the microstructural evolution and the deformation of simulation cells.

A. Simulation setup
The simulations were performed using cells comprising 1 024 000 atoms. Within each cell, a total of 80 × 80 × 80 bcc unit cells were present, with each unit cell containing 2 atoms. The interatomic potential developed by Ackland et al. [39] was used to describe the potential energy of the iron atoms.
Our objective was to simulate the condition of a material experiencing external uniaxial stress. However, in a simulation cell with periodic boundary conditions, there exists no free surface onto which external forces can be directly applied. To address this limitation, the simulation cells were deformed to simulate such a condition, thereby driving the Virial stress [40] to reach specific values. This approach, in principle, achieves equivalence to the application of external uniaxial stress, provided that the Virial stress is equivalent to the Cauchy stress.
Under uniaxial stress, the simulation cells were gradually deformed by initiating an energy minimization process using the conjugate gradient method. The orthogonal nature of the cell vectors was maintained throughout the deformation process. The cells were deformed until the stress in the z direction reached values of 0, ±0.00001, ±0.00002, ±0.00005, ±0.0001, ±0.0002, ±0.0005, ±0.001, ±0.01, ±0.02, ±0.05, ±0.1, ±0.2, ±0.5, or ±1 GPa, while the stresses in the x and y directions remained at zero. A positive stress value indicates tensile stress in the z direction, while a negative value indicates compressive stress.
All simulations were performed using LAMMPS [41]. Irradiation damage was simulated according to the CRA. Throughout the entire simulation process, simulation cell vectors are allowed to deform under the constraints of keeping the initial uniaxial stress unchanged, and of keeping the cell vectors orthogonal. For each uniaxial stress value, we performed three identical simulations using different random numbers as seeds for the stochastic process of choosing random atoms and moving them to random locations. Data are presented as the average of three samples.

B. Creation-relaxation algorithm
The meaning of DPA here is different from the conventional definition of NRT-DPA by Norgett, Robinson, and Torrens [5]. In CRA, the meaning of DPA is taken literally as it incrementally introduces damage in the form of Frenkel pairs by displacing random atoms to random positions. In some previous works, it was denoted as canonical DPA (cDPA).
Our results obtained by CRA are essentially simulating conditions where thermal diffusion of defects is suppressed, for instance, because the temperature is low or by the presence of impurities. Albeit the difference of definition, Derlet and Dudarev [29] and Warwick et al. [33,35] showed in Fe, W, and Zr, respectively, that the saturation behaviors of swelling due to accumulation of defects can be obtained quantitatively through rescaling, although the scaling relations are not straightforward. Quantities obtained from CRA simulations and experiments can also differ by a factor of 10 [31]. On the other hand, collision cascade simulations with low recoil energies (∼100 eV) can reach defect contents comparable to CRA [28], hence part of the mismatch can be attributed to the lack of defect recombination arising from higher-energy recoils present in realistic irradiation conditions. The CRA allows us to infer the principal mechanisms underlying certain phenomena at a fraction of the cost of conventional cascade simulations.
In the work by Derlet and Dudarev [29], CRA was implemented such that upon displacement of an atom, a dose of 1/N DPA is accumulated, where N is the number of atoms in the simulation cell. Frenkel pairs were created one by one, and after each creation of a Frenkel pair the atomic configuration was relaxed through energy minimization. Subsequent work using CRA [31][32][33][34][35] showed that the creation of multiple Frenkel pairs per relaxation step, corresponding to a dose increment of 10 −3 DPA per relaxation step or more, can speed up the simulation and have little effect on the results.
In our implementation, Frenkel pairs were generated to correspond with a dose increment of 10 −3 DPA per relaxation step. During each relaxation step, the atomic configuration was initially relaxed without altering the simulation cell's shape. Subsequently, the atomic configuration and the simulation cell were relaxed together until the Virial stress achieved the desired target values. This approach ensured numerical stability in the energy minimization procedure. The described procedure was repeated until all samples reached a cumulative dose of 3 DPA.

C. Defect relaxation volume density tensor
The concept of eigenstrain * was introduced by Mura [38]. Eigenstrain is a generic term referring to nonelastic internal strain. Reali et al. [37] proved analytically that the defect relaxation volume density tensor ω is equivalent to the eigenstrain of defects. In the absence of other sources of nonelastic strain, one can write * The defect relaxation volume density tensor can be defined as [37] where a i j is the relaxation volume tensor of a defect situated at R a . The relaxation volume tensor, together with the elastic constants, can describe the elastic properties of a defect in the asymptotic limit far from the defect. The values of i j of various defects can be obtained using density function theory calculations [42][43][44][45][46] and molecular statics [47,48].
Knowing the spatially varying ω(x), one can calculate the effective body force [37,49]: where C i jkl is the elastic stiffness tensor. This body force can be supplied to a continuum model, such as FEM, for calculations in the component scale.
Notably, component-scale simulations using FEM describe spatial scales many orders of magnitude larger than atomistic simulations. Assuming a component having a dimension of 1 m 3 , if one considers 1 million elements, the linear dimension of an element is 1 cm, which is much larger than atomic scale simulations using Å or nm as units. As a result, in engineering simulations, it suffices to obtain homogenized information pertaining to each element, without requiring intricate knowledge of the underlying microstructure. However, from a scientific perspective, our interest remains in comprehending how macroscopic properties evolve as a consequence of microstructural changes.
Following Mura [38], the total strain is expressed as the sum of eigenstrain and elastic strain: Assuming linear elasticity, the stress is related to the elastic strain through Hooke's law: Since eigenstrain equals to the relaxation volume density tensor, one can write the stress: Taking the average over the whole simulation cell, the average stress isσ After rearranging terms, one can writē where S i jkl is the elastic compliance tensor, which is the inverse of the elastic stiffness tensor C i jkl [50]. In this work, C i jkl is calculated from the interatomic potential. The average strain i j can be obtained according to the change of the simulation cell vectors. The average stressσ kl is the applied stress.
Our applied stress only has one nonzero component in the z direction. Using the compliance tensor appropriate for bcc crystal symmetry, the average relaxation volume density tensor can be simplified tō 073604-3 Using the above expression, we can obtain the values ofω as a function of uniaxial external stress at different levels of damage or dose. Such homogenized information can be applied to continuum models directly due to the different length scales being considered. The total defect relaxation volume equals the total defect formation volume, as the total number of atoms is conserved. Therefore, the magnitude of volumetric swelling is readily obtained from the trace ofω.

D. Crystal growth
The crystal growth under irradiation can be understood by examining the changes in simulation cell vectors [27]. One can write the matrix of cell vectors as In the case of a perfect lattice, this can be separated into two parts: where P is the primitive cell vector matrix, and N represents the repetition of the primitive cell. For the bcc lattice, the primitive cell vector matrix is By inverting P, one can get N = ⎛ ⎜ ⎝ n xy + n xz n yy + n yz n zy + n zz n xx + n xz n yx + n yz n zx + n zz n xx + n xy n yx + n yy n zx + n zy where n i j = L i j /a represents the repetition of primitive cells in Cartesian coordinates.
In an orthogonal cell, we can write where n xx , n yy , n zz is the number of unit cells in the x, y, and z directions, respectively. In the case of a purely elastically deformed simulation cell, the primitive cell vector matrix can be written approximately according to linear elasticity as where e = Sσ. On the other hand, we can also obtain the elastically deformed primitive cell from numerical calculations, which include nonlinear contributions, where and a x , a y and a z are the mean lattice parameters of an elastically deformed orthogonal unit cell in Cartesian coordinates. The three lattice parameters can be obtained from the deformed perfect simulation cell. We used this in our calculations.
When a simulation cell is subjected to irradiation, the crystal structure evolves. Assuming one can still identify the dominant crystal structure, we can estimate the number of repeating cells or planes along x, y, or z direction from: where L is the matrix of deformed cell vectors. The change of N over the course of irradiation is attributed to the plastic deformation caused by the accumulation of defects.

A. Change of macroscopic quantities
The change of the average strain¯ and the average relaxation volume density tensorω as a function of irradiation dose in units of DPA will be examined at different levels of uniaxial stress. The average strain under various conditions can be calculated from the simulation cell vectors by referencing the perfect crystal lattice under stress-free conditions at 0 DPA. Figure 1 shows the¯ of simulation cells subjected to different stress values at 0 DPA. The corresponding linear elastic approximations are also included in the plot. It is observed that even at stress levels as high as ±1 GPa, linear elasticity offers an accurate description of the elastic strain. This validates the use of Eq. (9) for estimatingω.
The change of¯ due to the application of external stress as a function of DPA is shown in Fig. 2 Similarly, one can observe the change ofω in Fig. 3. All components increase linearly up to about 0.01 DPA, and nonlinearly to about 0.03 DPA, regardless of the applied stress. We note that all values now start from zero, meaning there are initially no defects in the system. The linear increment at low dose is due to the generation of isolated SIAs and vacancies [29]. We will further verify this point in the following sections.
Since the elastic field of a defect does not need to be isotropic, external stress σ ext i j can affect the energy of a system containing a defect, where the elastic energy of a defect subjected to an external stress field [45] can be written as E el = −σ ext i j i j , where i j is the defect relaxation volume tensor. This means that the preferred orientation of defect clusters is biased by external stress, which affects the orientation of subsequently formed lattice planes. By inspecting the change ofω, we expect lattice planes to form preferentially with respect to the direction and value of stress. Mason et al. [31] show that in the case of a thin film, if the x and y dimension are constrained, lattice planes will form perpendicular to the z direction, which is similar to the tensile stress condition. We suggest a similar phenomenon is happening here, with further investigation following. Figure 4 shows the average strain¯ against external uniaxial stress. One can observe the¯ 11 and¯ 22 decrease and¯ 33 increase almost linearly with respect to stress. After eliminating the linear elastic effect, one can observe similar behavior ofω in Fig. 5. The main difference is that all the slopes are less steep. At 0 DPA, it is correctly observed that allω = 0.
At low dose below 0.2 DPA, components ofω show an almost linear relation to the external stress. However, when the dose is larger than 0.2, and when the uniaxial stress is small, within ±0.1 GPa, fluctuations in data can be observed. There are two possible reasons. First, the fluctuations mean that simulation cells deform without preferential direction. Second, we cannot observe the bias properly due to insufficient simulations. Both reasons imply that external stress is not strong enough to significantly affect the change of microstructure of current simulations.
We note that the dislocation density peaks at 0.2 DPA, at which point the microstructure contains larger dislocation loops due to the coalescence of smaller loops, at which point the dislocation network starts forming (we will discuss the microstructure below). The formation of dislocation networks reduces the ability of the system to respond to external stress since the energy barrier of changing the morphology of a dislocation network is much higher than for individual small loops and other localized defects.
The spatial fluctuation of stress could be inspected through the von Mises stress (VMS). VMS is often used in engineering to predict if a material will yield or fracture. The VMS σ vms is related to the Cauchy stress tensor σ i j by The atomic VMS is defined according to the atomic Virial stress, which allows us to examine the spatial distribution of VMS. Figure 6 shows the spatially average and standard deviation of VMS. We can see there is little difference between large compressive stress, tensile stress, and stress-free conditions. On the other hand, at a dose larger than 0.2 DPA, theω show clear preferential behavior only when the absolute value of external stress is larger than 0.5 GPa, at which point its magnitude is comparable to the VMS. One possible explanation is that in order to affect the microstructural evolution of a dislocation network as it is forming, external stress needs to be able to induce an energy change comparable to the energy scale of the network as a whole, rather than to the energy scale of individual defects. Consequently, the response of the system changes for a dose larger than 0.2 DPA.
We attempt to mapω to simple functions. This may help develop other continuum models, such as FEM. We assumē ω 11 andω 22 are degenerate, i.e.,ω 11 =ω 22 . We found thatω 11 andω 33 can be fitted to simple functions of dose and external stress, such that, where φ is the dose in units of DPA. After performing leastsquares fitting to all data, we get We acknowledge that Eq. (19) provides a simplified representation ofω solely depending on the dose and a single stress component. Generally, the behavior ofω is governed by multiple parameters, such as temperature, dose rate, and all components of the stress tensor σ. However, in this particular case, we are assuming irradiation conditions consistent with the athermal regime of microstructural evolution, where the dose rate is much higher than the rate for defect migration. Consequently, the temperature and time parameters are irrelevant [28,29].

B. Irradiation-induced swelling and creep
In general, swelling is represented by the volumetric strain where V and V 0 are current and initial volumes. Here the initial volume V 0 is that of the cell at 0 DPA deformed by external stress, not that of the stress-free cell. We plotted the change of the average atomic volume and the volumetric strain as a function of dose in Figs. 8 and 9, respectively. The swelling saturates to about 3% at a high dose, showing only minor dependence on stress, which is compatible with previous work by Derlet and Dudarev using other interatomic potentials [29].
One can see that even before irradiation, the atomic volumes are not the same at different stress values. This can be understood by considering that even in isotropic materials, the elastic deformation of the volume due to an externally applied uniaxial stress is in general not zero, because where E is Young's modulus and ν is Poisson's ratio. We observe that for ν < 0.5, when σ app is positive, the volumetric elastic strain represented by e v should be positive. This is consistent with what we can see in Fig. 8.
In fact, in the general case where a total strain with elastic and plastic components is acting on the cell, we can rewrite

Eq. 21 using (anisotropic) linear elasticity theory into
where V p is the volume of stress-free unirradiated perfect lattice cell. One can simplify it into where ω v = Tr(ω) =ω 11 +ω 22 +ω 33 is the volumetric eigenstrain. The elastic strain does not contribute to the volume change to first order, and higher-order terms are negligible in the linear elastic regime of deformation. Consequently, the volumetric eigenstrain is given by the plastic deformation only. See Fig. 10 for a plot of the eigenstrain, which at high dose is also found to saturate to about 3%. As expected, we can see that volumetric strain and eigenstrain are very close to each other.
Both v and ω v are mildly stress dependent. We inspected their dependence on stress in Figs. 11 and 12, noting a seemingly linear relationship where the swelling increases as tensile stress increases. Linear fits were performed to determine the slopes. We should note that ω i j is calculated assuming linear elasticity according to Eq. 9. The values of ω i j at ±1 GPa may not be very accurate due to a slight departure from the linear elastic regime, see Fig. 1. In Fig. 13, the slopes of linear fits of Fig. 11 and 12 and other similar data are shown. We can see the slopes saturate at about 8 × 10 −4 GPa −1 . The effect of stress on volumetric swelling exists but is small. From the experimental perspective, irradiation-induced creep is commonly separated into two parts; these are creep in the absence of swelling, and swelling-enhanced creep [22,24]. In the linear regime, these two contributions are described by a phenomenological equation: where p is the deviatoric plastic strain, σ app is the applied stress, B 0 is the deviatoric creep compliance for irradiation creep deformation, and D is the creep-swelling coupling coefficient for irradiation creep deformation. This equation applies to cases where the plastic strain increases as the dose increases, and swelling is comparatively small. According to Zinkle [22], experiments typically determine values for B 0 and D within the range of 0.2 to 0.45 of the melting temperature. The phenomenological equation is used to quantify creep under irradiation conditions where the thermal diffusion of defects is active. The expression is not an appropriate choice for describing athermal irradiation conditions, where thermal diffusion of defects is entirely suppressed. Here, we shall only present an order of magnitude estimate of the creep rate in the athermal limit in the linear regime, and compare the value of B 0 to values reported from experiments at finite temperatures. Since the value of D is related to the climb-controlled glide of dislocations involving the diffusion of point defects, it is reasonable to exclude its contribution within this context. First, we separate the eigenstrainω into two parts, namely the hydrostatic volumetric swelling ω v and volume conserved creep ω , such that which is shown in Fig. 14. The ω 33 is monotonic increasing or decreasing according to the sign of external stress. At high stress, it is clearly seen that ω 11 and ω 22 have opposing sign to ω 33 . The external stress is clearly responsible for creep, though we note that the plastic strains are not increasing or decreasing linearly; in the limit of high dose, they instead saturate. Only at a low dose below 0.2 DPA can we observe near-linear behavior. We plotted in Fig. 15 the value of ω 33 /σ 33 against the dose φ below 0.2 DPA. The slopes of the curves yield a creep rate of about 0.3 GPa −1 DPA −1 , compared to values of B 0 reported as 0.5 × 10 −3 GPa −1 DPA −1 in ferritic steels [3]. Our simulations in the athermal regime find a deviatoric plastic strain developing at a rate of three orders of magnitude faster than in conventional irradiation-induced creep experiments. Further, we note that the deviatoric plastic strain saturates at high dose, which is associated with the development of a steady-state microstructure saturated with immobile vacancies. In the steady-state, further irradiation does not generate more defects [28], and as such no more plastic deformation can occur in the absence of thermally activated mechanisms.
Regarding the numerical accuracy of the simulations, first we note that thermal diffusion of defects is entirely suppressed in the athermal regime of microstructural evolution: Defects with high migration barrier O(1 eV) (such as vacancies) are immobile, while defects with smaller migration barrier O(0.1 eV) (such as interstitials) can be driven to move by the stress fields developing in the system as defects accumulate [29,32,35]. The mechanisms underlying irradiation creep at low temperatures are therefore expected to be very different from those at intermediate or high temperatures. We expect our simulations to be qualitatively comparable to experiments in irradiation conditions at which vacancies are immobile, with a quantitative accuracy on the order-of-magnitude scale, as seen in similar work on tungsten [28,31]. Further, we note that we are simulating pure iron, and therefore alloying and impurity effects are omitted, which act to slow down defect mobility. We conclude that current results remain a plausible way to examine the underlying mechanism of swelling and  creep due to irradiation in the athermal limit, and should provide sensible predictions at temperatures below the onset of vacancy migration.

C. Microstructural evolution
The change of the microstructure under different external uniaxial stresses was examined. The analysis and visualization of the samples were performed using OVITO [51].
First, the content of point defects is analyzed. Point defects were identified by Wigner-Seitz analysis, which is implemented in OVITO lattice site, one can count the number of atoms corresponding to the displaced configuration. The reference configuration is the initial defect-free simulation cell at 0 DPA. A Voronoi volume having no atom is considered as a vacancy, while a volume with two atoms occupancy is considered as an SIA. Voronoi volumes with three atoms occupancy or higher are rare, so we can safely discard them.
Similar to the work by Derlet and Dudarev [29], we further identify the isolated and nonisolated SIAs, but we used different criteria. We classify an SIA as belonging to a cluster if there is another SIA within a specified cutoff distance. We used a cutoff of 3.2 Å, which is a distance between the second and third nearest meighbors of an atom in the perfect iron bcc lattice. Since Wigner-Seitz analysis places SIAs on reference lattice sites, which are lattice points of the perfect crystal, our results remain the same for any cutoff distance between a and √ 2a, or 2.9 Å and 4.0 Å, respectively. In the work of Derlet and Dudarev [29], they identified the isolated SIAs if there is no other atom sitting within a cutoff of 0.93 of √ 3a/2, which is approximately 2.3 Å for bcc iron. However, we should remember that two atoms sit within a Voronoi volume of a reference lattice site, it only has 1 SIA. It is ambiguous to tell which one is the SIA, and so is the atomic position. If both are considered as SIA, it is double counting the total number of SIA. Instead, if we only consider the reference lattice site, we know the exact location of an SIA. Figure 16 shows the vacancy and isolated SIA contents. Interestingly, we can see there is little difference between samples at ±1 GPa and stress-free conditions. All of them reach a maximum isolated SIA concentration at around 0.03 DPA. The vacancy and isolated SIA concentrations start diverging at around 0.01 DPA, where larger defect clusters such as dislocation loops and SIA clusters are forming. At high doses, the vacancy and isolated interstitial contents reach a steady state with concentrations of 4.2% and 0.8%, respectively. They are in quantitative agreement with work by Derlet and Dudarev [29] using two other interatomic potentials.
Then, we analyzed the Laves phases cluster size. We identified the Laves phases using the polyhedral template matching (PTM) method [52]. According to Chartier and Marinica [36], one could identify the C15 Laves phases by finding the icosahedral structure. Essentially, C14, C15 and C36 Laves phases are made of Z12 and Z16 Frank-Kasper clusters. Z12 means icosahedral arrangement with 12 neighboring atoms. Z16 is a similar arrangement but with 16 neighboring atoms. PTM can only identify icosahedral arrangements, but since C15 is expected to stabilize in irradiated iron, they argued that the identification of Z12 is equivalent to C15, where independent analysis [53] confirmed this. We note that another method to identify different Laves phases was developed by Xie et al. [54]. Cluster analysis is used for the analysis of the sizes of Laves phases clusters. Cluster analysis uses a cutoff distance similar to the coordination analysis. For any given particle, they are iteratively added to a particular cluster. Finally, we got the cluster sizes and know which particle belongs to which cluster.
In Fig. 17, we plotted the cluster size against the dose in units of DPA. The counting of clusters is shown as color 073604-13 map. One can observe a bright spot at around 0.1 DPA, corresponding to the peak quantity. Our results are in qualitative agreement with the work by Chartier and Marinica [36]. We obtained a larger number of Laves phases clusters, where the sizes of clusters are smaller than the results by Chartier and Marinica. A possible reason is that they are using a different interatomic potential and working at finite temperatures. Figure 18 shows the volume concentration of clusters as a function of dose. Similarly, we see a peak in cluster concentration at around 0.1 DPA, after which the concentrations drop and eventually reach steady values. It is because some of the clusters are transforming into dislocation loops [36]. Interestingly, we found external stress to have little effect on the cluster content.
Finally, we analyzed the dislocation density. Dislocation analysis [55] being implemented in OVITO is used to identify the length and nature of dislocations. Fig. 19 shows the dislocation density for different Burgers vectors. The dislocation content peaks at around 0.2 to 0.25 DPA, after which it drops gradually. Only dislocations of a 2 111 and a 100 type are observed, where a 2 111 -type dislocations have significantly higher content. These observations agree with similar simulations [29,36]. The relative formation energy of a 2 111 and a 100 dislocations and C15 clusters [56] explains the relative densities as well as the point at which Laves phases cluster concentration starts to drop while dislocation density starts to rise, where Laves phases clusters nucleate dislocation loops [36,56].
Again, we do not see any evidence that external uniaxial stresses in the range of ±1 GPa affect the dislocation density and type in comparison to stress-free conditions. Figure 20 shows the microstructural content of isolated SIAs, vacancies, dislocations and Laves phases clusters sideby-side. It gives us a better idea of the spatial distribution of defects. Since we found they behave similarly under different external stress, we only plotted the stress-free samples. We can observe clearly that the content of isolated SIAs drops and vacancies increases starting from 0.01 DPA. On the other hand, the Laves phases clusters peak at 0.1 DPA, where dislocations peak at 0.2 DPA. It is corroborated by the explanation that Laves phases clusters are transforming into dislocation loops [36].
According to what we found here, it appears that the microstructure in relation to isolated SIAs, vacancies, Laves phases clusters, and dislocations is not much affected by the external uniaxial stress. These defects are not the main cause of the change and polarization of theω under external uniaxial stresses.

D. Anisotropic planes formation
The anisotropy ofω is indicative of the polarization of defects formed during irradiation. As we found the content of isolated defects to be unaffected by the application of external uniaxial stress, we expect the anisotropy ofω to be caused by anisotropic crystal growth. Theω of different stresses start diverging at around 0.2 DPA, at which point dislocations start coalescing, forming dislocation networks and eventually lattice planes [29]. Experimentally [24], the number and size of dislocation loops in AISI 316 stainless steel under stress have been studied. The loop number density increases in the lattice planes normal to the applied stress, with a decrease in the number of loop orientations parallel to stress, but the loop size distributions are not changed by applied stress. This suggests microstructural evolution under irradiation can be affected by external stress. In our case, we argue that the external applied stress leads to a preferential orientation of lattice plane growth, leading to the anisotropic changes inω.
To test this hypothesis, nonisolated SIAs were identified and visualized, as shown in Fig. 21. It is evident that under tensile stress, planes form perpendicular to the direction of uniaxial loading. This observation aligns with the results 073604-15 obtained from CRA simulations conducted by Mason et al. [31,32], where the formation of interstitial planes normal to the z direction was observed when the cell dimensions were constrained in the x and y directions to mimic a scenario where the sample is attached to a fixed substrate. Conversely, when our samples were subjected to compressive stress, the orientation of plane growth exhibited less obvious patterns.
The quantification of changes can be achieved by considering the relationship between cell vectors and the repetition of unit cells. Figure 22 illustrates the alterations observed in n xx , n yy , and n zz . Under a tensile stress of 1 GPa, the value of n zz increases from 80 to slightly above 83. This can be compared with the findings depicted in Fig. 21, where approximately three planes are observed. In the case of compressive stress, both n xx and n yy reach 82, while n zz decreases to 78.5. It can be inferred that planes are forming in a manner that contributes to growth in the x and y directions while experiencing a reduction in the z direction.
Based on the data presented, it can be concluded that the anisotropy ofω under external stress can be attributed to the stress-dependent orientation of the interstitial planes that are formed.

IV. CONCLUSION
The influence of external uniaxial stress on athermal irradiated iron was investigated using atomic-scale simulations employing the creation-relaxation algorithm. This approach enabled the examination of irradiation-induced swelling and creep by inspecting the changes in defect relaxation volume density (or eigenstrain) at high irradiation doses.
Analysis of the eigenstrain of defects revealed anisotropic changes when external stress is applied. The polarization effect was observed to commence at approximately 0.1 to 0.2 DPA, where self-interstitials form dislocation loops, followed by the subsequent development of dislocation networks.
Further investigation of the microstructure unveiled that the anisotropic changes in the eigenstrain were attributed to the anisotropic growth of interstitial planes driven by the external stress. The application of external stress exhibited minimal impact on the density of isolated SIAs, vacancies, Laves phase clusters, and dislocations, and vice versa.
However, under tensile stress, a preferential formation of planes with normals parallel to the loading direction was observed, while compressive stress resulted in a more complex pattern of planes with normals perpendicular to the loading direction.
The obtained findings elucidated the significant role of external stress in the behavior of heavily irradiated iron. The simulations provided valuable insights into the underlying mechanisms governing irradiation-induced swelling and creep in iron. These insights hold the potential to inform the design and optimization of advanced reactors.