Dislocation loop bias and void swelling in irradiated α -iron from mesoscale and atomistic simulations

Dislocation loops are ubiquitous in irradiated materials, and dislocation loop bias plays a critical role in void swelling. However, due to complicated interactions between dislocation loops and point defects, it is challenging to evaluate the bias factors of dislocation loops. Here, we determine the bias of sessile < 100 > loops in α -iron using a recently developed atomistic approach based on the lifetime of point defects. We establish a mechanistic understanding of the loop interaction based on the diffusion tendency of point defects near the loop core region. Mobile self-interstitial atoms tend to be absorbed from the edge of the loop, and a trapping region perpendicular to the habit plane of the loop exists. The dislocation loop bias is found to be substantially lower than those of straight dislocations in α -iron and should be included in swelling rate estimates. With the obtained sink strength and bias values, agreement is achieved with experimental results for both absolute values and temperature dependence.

T he unit interactions of defects at atomistic scales during irradiation processes [1][2][3][4] could significantly affect macroscopic material properties, such as void swelling [5][6][7][8][9][10][11] , degrading mechanical performance, and limiting the lifetime of structural components in fission and fusion reactors 12,13 .The generally accepted mechanism for void swelling is the preferential absorption of mobile interstitials over vacancies at certain sinks, defined as bias 12,14,15 .The bias factor values for different types of sinks are essential parameters for modeling void swelling and overall microstructural evolution.However, it is fundamentally challenging to study the underlying complex defect interaction processes between sinks and point defects, which are too fast to observe using experimental tools like transmission electron microscopy (TEM) and too slow to model by conventional atomistic simulations [16][17][18][19][20] .Moreover, some previous theoretical models, e.g., in α-iron, have overestimated bias values by an order of magnitude compared with those estimated from experimental results 10,21,22 .Therefore, to precisely evaluate the bias in structural materials, a fundamental understanding of interactions between sinks and radiation-induced defects is essential.
One of the dominant microstructural features of sinks for mobile point defects in irradiated materials is dislocation.Dislocation bias factors have been widely studied using elasticity theory 9,10 , production bias model 23,24 , dipole tensor method 22,25 , and numerical approaches [26][27][28][29] .However, it is still challenging to achieve close agreement with experimental results of swelling rates based only on dislocation bias.
Dislocation loops are also ubiquitous in irradiated α-iron [30][31][32][33][34][35][36][37][38] and are known as biased sinks for mobile point defects 39 , including interstitials and vacancies.Comparatively, dislocation loop bias is much less frequently studied due to the extremely complex three-dimensional strain field near it 40,41 .For instance, the formalism of the strain field around a dislocation loop is daunting even in an infinite isotropic material 42,43 , and it is challenging to apply these complex strain field equations in the bias calculations.Therefore, different approximations were made in previous theoretical efforts.Brailsford et al. 15,44 originally developed a rate-theory framework for describing void swelling and assessed the role of dislocation loop bias.This model was employed by Bullough and Woo et al. 39,41 to calculate sink strength and bias factors of dislocation loops with the effective medium approximation.In these studies, the dislocation loop is considered a spherical sink.Dubinko et al. 45 employed a toroidal internal boundary condition for dislocation loops, which emphasized the importance of sink topology.The description of dislocation loop topology was improved by Rouchette et al. 40,46 through the phase-field method 47 , which determined sink strength more accurately by considering the anisotropic elastic interactions between point defects and loops.Nevertheless, there have been two significant limitations.First, the effects of the core region near the dislocation loop cannot be accurately modeled because the atomistic details of loops were not considered.Second, the influences of dislocation loops on defect transport properties, e.g., migration energy barriers (MEBs), are inaccurate since the interactions between point defects and dislocation loops are stronger than those between point defects and straight dislocations, due to the larger strain field exhibited by the loops 46 .This will result in regions where the actual saddle points differ noticeably from the estimations based on the isotropic/anisotropic elasticity theory approaches, subsequently impacting calculations of capture efficiencies, sink strengths, and loop bias.
In this study, the dislocation loop bias is calculated based on the lifetime of point defects.Dislocation loops with radii from 1 nm to 5 nm are studied based on the experimentally reported loop sizes in neutron-irradiated iron 31 .Only sessile < 100 > loops are considered since they are the dominant type in bcc iron when the temperature is higher than 548 K 31,48,49 .Spontaneous absorption region (SAR) is calculated based on the binding energy using Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) 50 , in which the mobile point defects are considered absorbed by the loop.The MEBs of point defects at every lattice point near the dislocation loop are calculated using   71 is shown in pink.Note that the MEBs of saddle points (ii and v, i and vi, iii and viii, iv and vii in Fig. 1e) look identical, although there exists a very small numerical difference.
the self-evolving atomistic kinetic Monte Carlo (SEAKMC) package [51][52][53] .The statistical diffusion tendency considering all different dumbbell configurations is analyzed at each atomic position close to the loop.The point defect lifetime is then calculated using atomistic kinetic Monte Carlo (AKMC), with which the capture efficiency, sink strength, dislocation loop bias, and swelling rates are also calculated.Temperature, loop size, and loop density effects are systematically studied.The obtained bias factors are compared with previous theoretical results of loop bias and straight dislocation bias.Estimated swelling rates and experimental data are also compared, and the effects of various factors that influence void swelling in irradiated α-iron are discussed.

Results and discussion
Migration energy barriers (MEBs).Figure 1 shows the MEBs calculated based on obtained saddle point configurations (illustrated in Fig. 1e with blue spheres labeled from i to viii, respectively) associated with the first nearest neighbor diffusion of the [110] dumbbell.For instance, the MEBs toward ½ 1 11 direction through saddle point ii are given in Fig. 1a.It is found that MEBs are highly anisotropic and exhibit a complex dependency on the position relative to the dislocation loop due to the previously mentioned complicated strain field of the dislocation loop.For the [110] dumbbell, the region with lower MEBs is primarily along the normal direction of the loop habit plane, and the region with higher MEBs circles the loop.In addition to the change in energy barriers, note that dumbbell diffusion via a given saddle point leads to a different final configuration from the initial configuration.Taking a [110] dumbbell as an example, the final configuration will turn into a ½10 1 dumbbell through saddle point ii or v.Besides the [110] dumbbell, all five other dumbbell configurations present distinct MEB dependencies, which have been shown in Supplementary Fig. 1 together with the corresponding final configurations through different saddle points.Since the strain field is not homogeneous and since previous researchers 39,45 also stated that the loop bias was found to depend on the loop nature and the point defect shape at the saddle point, the spatial dependence of MEBs in the vicinity of the dislocation loop is thus fully considered, which is challenging using the conventional rate theory approaches.
We find the cumulative effects of MEBs and the configurations of dumbbells contribute significantly to the bias factor.We conduct comparative simulations to demonstrate these effects (see data in Table 1, which shows the bias factor (0.0089) of a 5 nm loop at 573 K at a dislocation loop density equal to 10 9 cm −2 ).For comparison, a significant overestimation (6.8 times higher) is observed if the MEBs calculated using SEAKMC are replaced by those calculated using only the dipole tensor method.Correspondingly, we notice an obvious underestimation of capture efficiencies for both self-interstitial atoms (SIAs) and vacancies using the elasticity method.This significant difference is caused by the cumulative effect of the huge number of MEBs and the number of diffusive steps of each trajectory (e.g., about 0.7 million diffusive steps when the dislocation loop density equals 7 * 10 9 cm −2 at 573 K), although the differences in MEBs at most of the atomic positions are small (details in Supplementary Fig. 2).We further show the effect of dumbbell configurations by considering only one of the six dumbbell configurations in the model.We find two different bias factors for different dumbbell configurations due to crystallographic symmetry.Using either ½110 or ½1 10 configuration leads to a negative bias factor around −0.01; ½101, ½10 1, ½011, and ½01 1 dumbbells result in a much higher value, which is about an order of magnitude higher than the value from our model.Thus, it is essential to include all possible initial configurations and the configuration change at each diffusive step in KMC simulations.
Point defect transport mechanism near dislocation loop.Near the dislocation loop core region, SIA diffusion shows two pathways, illustrated in Fig. 2a.In Path I, a mobile SIA near the loop habit plane tends to be absorbed directly by the edge of the loop.In Path II, a SIA far from the loop habit plane migrates toward the habit plane region along the curved arrow paths.This transport behavior in Path II results from a strong repulsive interaction between the loop and the SIA.Comparatively, for vacancy transport, only one path is found, wherein a mobile vacancy tends to be absorbed by the loop from ~35°to the habit plane, as shown in Fig. 2b.
This point defect transport mechanism is revealed by analyzing the expected diffusion direction of a point defect at every atomic position, termed as diffusion tendency (details in Supplementary Note 2). Figure 2c shows the diffusion tendency of a SIA surrounding a 2 nm [100] loop, where the SARs are labeled with grey spheres.The (110) plane cross-section is shown since it has the highest atomic density, and the SIAs on this plane primarily show diffusion tendencies parallel to this plane.Four regions with different SIA transport behaviors are observed and labeled by α, β, γ, and δ.Region α is an attractive region (red arrows) indicating that a [100] interstitial loop generally absorbs SIAs from the edge of the loop, which corresponds to Path I. Region γ represents a repulsive region (blue arrows), and SIAs in this region tend to diffuse away from the loop.The fan-shaped region located near region γ is region β.Although the region β is attractive, SIAs in this region are prevented from being absorbed by the loop directly due to the strong repulsion in region γ.This phenomenon leads to the formation of a special region, labeled δ.Since most of the SIAs tend to diffuse towards the center of δ, they can be regarded as being trapped in δ and have a high chance to stay there for extra diffusive steps.This causes an increase in the lifetime of SIAs (τ i ) and lower absorption efficiency of SIAs.The trapping phenomenon in δ, the attractive influence in β, and the repulsive influence in γ lead to curved paths for the diffusion of SIAs, which corresponds to Path II.In comparison, Fig. 2d shows the diffusion tendency of a vacancy, where the attractive interaction between the loop core and vacancies is ~35°to the habit plane (red arrows).Vacancies are easily absorbed by the loop center in region γ, while they are repulsed in region β.Recent research by McElfresh et al. 54,55 using a phase field method has shown a similar behavior of vacancies approaching interstitial loops in Mo, which is that vacancies diffuse toward the loop from compressive stress gradients rather than directly from the habit plane, demonstrating highly anisotropic diffusion of point defects in response to the spatially complex stress fields created by the dislocation loops.The observations of point defectloop interactions are commensurate with the calculated stress/strain field of the loop (details in Supplementary Fig. 3).This mechanism is also supported by the lifetime of SIAs (Fig. 2e) and vacancies (Fig. 2f) around the loop.In Fig. 2e, the lifetime of a SIA in α is comparatively shorter than that of other regions with a similar distance to the loop surface.For a SIA in γ, the lifetime pattern indicates the diffusion trajectories are much longer than those being directly absorbed by the loop, despite the proximity of region γ to the loop center.In Fig. 2f, the lifetime of vacancies at ~35°to the habit plane is much shorter than those along other directions, which agrees with the diffusion tendency analysis.Capture efficiency and bias of < 100 > loops.Capture efficiency and bias of loops with radii from 1 nm to 5 nm are calculated with the obtained lifetimes.We find that loops with radii between 2 nm and 5 nm show lower bias factors (≤ 0.02 at 573 K, details in Supplementary Fig. 4) than that of a 1 nm loop (0.17 at 573 K).Since loops with sizes larger than 5 nm are commonly observed in irradiated iron 31,34,35 , the results of a 5 nm [100] loop as a function of dislocation loop density are shown in Fig. 3 as an example.The capture efficiency of a SIA (Z i , Fig. 3a), a vacancy (Z v , Fig. 3b), and the bias factors (Fig. 3c) decrease with decreasing loop density, which is consistent with previous calculations [39][40][41]45,46 . At a iven temperature and loop density, Z i is greater than Z v , due to the respective difference in the interaction energy 56 .The values of both Z v and Z v decrease with increasing temperature from 573 K to 873 K; however, the bias factors increase with rising temperatures.This trend differs from the calculation of vacancy loop bias by Woo et al. 39 and dislocation bias by Chang et al. 26,27 , Bakaev et al. 28 , and Seif et al. 22 .This phenomenon is attributed to the trapping phenomenon discussed in Point defect transport mechanism near dislocation loop.Since the kinetic energy of the system increases with temperature, it is easier for SIAs to overcome the MEBs, reducing the time a SIA remains in the trapping region (i.e., reducing the lifetime τ i ).In contrast, there is no trapping region for vacancies.
The obtained dislocation loop bias is substantially lower than that of straight dislocations calculated by previous studies 10,22,[26][27][28][29]57,58 , as shown in Fig. 4, even screw dislocation bias calculated using the same approach 29 . Note hat the bias value of a screw dislocation from Bakaev et al. 28 is close to the loop bias in this study.However, this may be a coincidence since Chang et al. 26 and Bakaev et al. 28 also reported negative bias for screw dislocations under other conditions.For example, these studies have reported bias values between −0.04 and −0.06 in the atomistic-analytical screw dislocation case without external stress, which may indicate systematic underestimation.Bias factors of loops in this study are also significantly lower than those of edge dislocations, either compared with calculations based on the driftdiffusion model and the elasticity theory by Wolfer 10 , Seif et al. 22 , and Borodin et al. 57 , or those calculated using atomistic approaches provided by Bakaev et al. 28 , Chang et al. 26 , and Kohnert et al. 58 .
Comparisons with experimental results.The estimated swelling rates are compared with experimental results [31][32][33][34][35][36][37][38]59,60 , shown in Fig. 5, which accounts for the sink strength of dislocation loops, straight dislocations (screw dislocations, data from Hao et al. 29 ), and voids. The dnsity of loops, dislocations, and voids, as well as the void sizes in the calculation, are based on the experimental results of Horton's work 31 since it provides very detailed data of microstructures.Both the calculations of this study and experimental results show that the void swelling rate increases with increasing temperature until reaching a peak value (at ~673 K) and declines with higher temperatures after the peak 31,32 .Before reaching the temperature where peak swelling occurs, small loops (5 nm-21.5 nm in radii) have been observed in ref. 31 as the dominant microstructure feature with high densities (> 10 20 m −3 ), indicating that dislocation loops play a critical role in void swelling within this temperature range.The calculated swelling rates increase from 0.008 % dpa −1 to 0.049 % dpa −1 , which is consistent with experimental data from Horton et al. 31 and Budylkin et al. 33 , corroborating the increasing loop bias.Alternatively, the contribution of straight dislocations to the void swelling changes little in this temperature range since the dislocation density (~10 14 m −2 ) and structure change little before reaching the temperature where peak swelling occurs 31 .
At 673 K, the calculated swelling rate from this study is ~0.044 % dpa −1 , which is consistent with the average value of several experimental studies [31][32][33][34][35] , although the residual impurities 32 , irradiation dose rates 33 , and cold-work level of samples 34,35 may cause fluctuation of results in individual experimental studies.When the temperature exceeds ~700 K, the vacancy emission dominates the void swelling, significantly limiting the void growth in experiments 61 .Based on our calculations (details in Supplementary Note 6), voids will lose 75.1% of their volumes at 723 K, because of which the bias of sinks is no longer the leading factor for the overall swelling rate, although loop bias factors increase with increasing temperature.The calculated temperature dependence over 673 K is also consistent with experimental results [31][32][33]36 , which is a significant improvement compared with previous theoretical studies.

Conclusions
The loop core effects and the cumulative effects of millions of MEBs significantly influence the defect transport properties of point defects.Based on the obtained MEBs and diffusion tendency, a mechanistic understanding of point defect diffusion near a dislocation loop in α-iron is developed.We identify two diffusion paths of a SIA near the dislocation loop, depending on the initial position of the defect.For Path I, a mobile SIA near the habit plane of the loop tends to be absorbed directly by the edge of the loop.For Path II, a SIA away from the habit plane of the loop diffuses along curved pathways toward the habit plane region.This mechanism results in trapping regions near the loop, which prolongs the lifetime of a SIA and leads to lower capture efficiency, which subsequently impacts the sink strengths and bias values of dislocation loops.
The < 100 > interstitial loop bias factors are lower than 0.1 with radii between 2 nm and 5 nm, substantially lower than those of straight dislocations.This conclusively demonstrates the significance of including atomistic details of the interaction between point defects and loops to void swelling.The calculated swelling rates increase with rising temperatures before peak swelling and decrease after 673 K due to the domination of vacancy emission from voids.These results confirm we have considered the decisive aspects of void swelling in irradiated α-iron across varying temperatures.
Finally, the exceptional agreement with experimental results constitutes extraordinary progress in swelling performance prediction and microstructural evolution simulation in irradiated materials.The generality of this approach allows straightforward extensions to provide insights into a variety of other situations as well, for example, modeling the interaction of microstructures and point defects in other structural materials, such as fcc metals or alloys.

Methodology
The lifetime, capture efficiency, and bias factor of point defects.The lifetime of a defect is defined as the time required for the point defect to diffuse from an initial position to a sink, which can be a dislocation loop, a straight dislocation, a void, etc.With the lifetime of a point defect, the capture efficiency can be defined following the approach derived in Hao et al. 29 , where τ α represents the lifetime of the point defect α, e.g., a SIA or a vacancy, under the influence of a sink, denoted as real walk in this paper.In comparison, the superscript r refers to random walk diffusion, which is the lifetime of a point defect without the interaction with the sink.The capture efficiency is the ratio of τ r α and τ α , which quantifies the preferential diffusion of a point defect to a specific sink with reference to bulk diffusion.The lifetime of the point defect is found to be proportional to the volume of simulation system, as demonstrated in Supplementary Fig. 5.This relationship can be utilized to predict the point defect lifetime in larger systems.The tendency of efficient absorption by the sink corresponds to a higher value of Z α .Following Golubov et al. 56 , the bias factor B of a specific sink is determined from the capture efficiencies of SIAs and vacancies using the following equation in which the subscripts i and v stand for SIA and vacancy, respectively.The bias factor quantifies the preferential absorption of SIAs compared with vacancies for a specific sink.Therefore, the bias factor can be calculated by substituting capture efficiencies with lifetimes: Swelling rate calculation.Swelling rate calculation is derived in refs. 21,56and calculated based on the sink strength for comparison with experimental results where ε is the survival fraction of point defects, which means the fraction of the Frenkel pairs escaping the in-cascade recombination, dependent on materials and irradiation conditions, usually between 0.1 and 0.3 [62][63][64] .This term provides the number of mobile point defects that should be considered in the diffusion calculations.The survival fraction ε is selected to be 0.1 in this study based on earlier atomistic simulations of displacement cascades and subsequent annealing in bcc iron 64,65 .k 2 represents the sink strength in the ideal system, subscript V stands for the void, l for the dislocation loop, and d for the straight dislocation.J em v is the vacancy emission rate 56,66 52,64 .The MEBs far from the loop core are calculated using the dipole tensor method 68,69 .As designed, the differences between MEBs on the edge of the box calculated using SEAKMC and the dipole tensor method are less than 0.002 eV.Detailed information is provided in Supplementary Note 1.
Atomistic kinetic Monte Carlo (AKMC).AKMC is employed to calculate the lifetime of a point defect.The setup of AKMC is shown in Fig. 6.A spherical system is set with the dislocation loop located at the center of the sphere, which is equivalent to uniform distribution of loops.The SAR is defined in this work as the region where the sink spontaneously absorbs the point defect.
Areas in red and grey are the regions in which the MEBs are calculated by SEAKMC and the dipole tensor method, respectively.Random boundary conditions 29 are employed to reduce the anisotropic impacts at low temperatures and high dislocation loop density simulations.Following this setup, the lifetime of a point defect is the statistical average diffusion time starting from the boundary to the SAR.The statistical accuracy is provided in Supplementary Fig. 6.
Loop construction and spontaneous absorption region (SAR) using LAMMPS.In bcc iron, < 100 > interstitial loops are the dominant type of dislocation loop in irradiated iron 31,70 .The configuration of a [100] interstitial loop is constructed by inserting two layers of SIAs into a perfect bcc system with the Burger's vector along the z-axis ([001]), while the x-axis is [100], and the y-axis is [010].The interatomic potential developed by Ackland et al. 67 is employed for structural relaxation using LAMMPS 50 .SARs have been determined based on the binding energies of point defects using the same potential.The criteria are set to be lower than −0.1 eV and −0.05 eV for an SIA and a vacancy, respectively.The definition of binding energy is the same as the interaction energy in Chang et al. 26 , which is the difference in the formation energy with and without the loop: The SAR for random walks of SIAs and vacancies is the same in this work as the reference state for real walks.

Fig. 1
Fig. 1 Migration energy barriers of the [110] dumbbell saddle points in bcc iron.a-d Each contain the MEB map of two saddle points given by blue spheres in (e).Figure (a) for the saddle point ii and v, (b) for i and vi, (c) for iii and viii, and (d) for iv and vii.Axes x, y, and z in (a) stand for < 100 > , < 010 > , and < 001 > directions.White spheres in (e) represent the [110] dumbbell.MEBs close to bulk values (± 0.001 eV) are intentionally hidden.The dislocation loop calculated by dislocation analysis (DXA) by Ovito71 is shown in pink.Note that the MEBs of saddle points (ii and v, i and vi, iii and viii, iv and vii in Fig.1e) look identical, although there exists a very small numerical difference.
Fig. 1 Migration energy barriers of the [110] dumbbell saddle points in bcc iron.a-d Each contain the MEB map of two saddle points given by blue spheres in (e).Figure (a) for the saddle point ii and v, (b) for i and vi, (c) for iii and viii, and (d) for iv and vii.Axes x, y, and z in (a) stand for < 100 > , < 010 > , and < 001 > directions.White spheres in (e) represent the [110] dumbbell.MEBs close to bulk values (± 0.001 eV) are intentionally hidden.The dislocation loop calculated by dislocation analysis (DXA) by Ovito71 is shown in pink.Note that the MEBs of saddle points (ii and v, i and vi, iii and viii, iv and vii in Fig.1e) look identical, although there exists a very small numerical difference.

Fig. 2
Fig. 2 The point defect transport mechanisms.Diffusion paths of SIAs (a) and vacancies (b) are both provided together with the repulsive regions near the loop.The diffusion tendency of SIAs (c) and vacancies (d) are calculated at each atomic position.The lifetime of SIAs (e) and vacancies (f) are calculated at 573 K and the loop density is 2.24 * 10 10 cm −2 (equals 1.78 * 10 22 m −3 when the loop radius is 2 nm).The SARs are labeled in grey atoms in the center of the simulation system for (c−f).

Fig. 3
Fig. 3 Calculated capture efficiency and bias as a function of the loop density.a Capture efficiencies of interstitials.b Capture efficiencies of vacancies.c Bias factors of a 5 nm [100] loop.

Fig. 4 A
Fig. 4 A comparison of loop bias in this study and straight dislocation bias from previous theoretical calculations.Data points close to 573 K are selected.The bias factors of screw dislocations are labeled by filled markers and edge dislocations by hollow markers.The (*) mark indicates a mixture of both dislocations, whereas screw dislocations are dominant (80% screw -20% edge).

Fig. 5 A
Fig.5A comparison of swelling rates in bcc iron between this work and experimental results.Filled and hollow symbols stand for neutron and ion irradiations, respectively.ρ l stands for loop density (m −3 ), ρ d for dislocation density (m −2 ).Microstructure information comes from Horton et al.31 .

Fig. 6
Fig.6The schematics of the AKMC simulation setup.The dipole tensor method is employed in the grey region, and the SEAKMC is employed in the pink region.SAR is shown as the blue circle in the center.Black spheres show an example of the point defect diffusion trajectory.

Table 1
Loop bias and capture efficiencies under different simulation conditions.
½110 dumbbell and ½1 10 dumbbell have almost the same results, the other four dumbbell configurations have similar results, which are different from those from ½110 dumbbell.The temperature is 573 K and the loop density is 10 9 cm −2 (equals 3.18 * 10 20 m −3 when the loop radius is 5 nm) for all these simulations.
67lk values equal to 0.3075 eV and 0.6401 eV for an SIA and a vacancy, respectively, based on the potential developed by Ackland et al.67.MEBs close to the loop core (about 58.54 million saddle points in a 23 * 23 * 23 nm 3 simulation box for a 5 nm loop) are calculated using SEAKMC , and N is the number of vacancies contained in a void.Migration energy barrier (MEB) calculation.Point defect migration to the first nearest neighbors is considered for MEB calculations, with