Cascade fragmentation: deviation from power law in primary radiation damage

ABSTRACT The sizes of defect clusters, produced in materials by energetic ion or neutron impacts, are critically important input for models describing microstructural evolution of irradiated materials. We propose a model for the distribution of sizes of vacancy and self-interstitial defect clusters formed by high-energy impacts in tungsten, and provide new data from in situ ion irradiation experiments to validate the model. The model predicts the statistics of sub-cascade splitting and the resulting distribution of primary defects extending over the entire range of cluster sizes, and is able to provide initial conditions for quantitative multi-scale simulations of microstructural evolution. GRAPHICAL ABSTRACT IMPACT STATEMENT We present a model, parameterized for tungsten, for the distribution of defect sizes in primary radiation damage, as an essential step in multi-scale modelling of microstructural evolution in irradiated materials.

Energetic neutron and ion irradiation changes mechanical and physical properties of materials, in part through the accumulation and evolution of primary defects, including vacancies, self-interstitial atoms (SIAs), and nano-sized dislocation loops, which form as a direct result of collision cascades initiated by the incident particles. Exposure to radiation is often measured using the notion of dpa (displacements per atom) or, equivalently, the number of point defects formed. However, this simple number provides no information on defect cluster sizes. In reality, the size statistics of defect clusters represents a critical factor determining how the microstructure of the material and its properties evolve under irradiation. The strength of the elastic interaction between the defect clusters as well as the lifetime of clusters, determined by the rate of evaporation and absorption of migrating point defects, depends entirely on the cluster size and separation. The initial distribution of defect sizes thus provides critical input in microstructural evolution models such as cluster dynamics or rate theory methods.
In tungsten (W), which due to its high thermal conductivity and low sputtering yield is the main candidate material for plasma facing components of tokamak fusion devices, including ITER [1], dense collision cascades produce large defect clusters [2]. These are readily visible in an electron microscope [3], making direct comparison between simulation and experiment possible. It has recently been discovered, both by direct molecular dynamics (MD) simulations [4] and transmission electron microscopy (TEM) observations [3], that the distribution of defects as a function of their size is well described by a power law where N is the size of a defect cluster (in terms of the number of point defects constituting the cluster), and S is a scaling exponent. A is a pre-factor proportional to the experimental or simulated defect cluster production rate. Events described by the tail end of the power law distribution are statistically significant, in contrast to, e.g. the tail end of a normal distribution. Hence a correct treatment of the upper size limit is as important as the main characteristics of the distribution itself. Physical considerations make it clear that the defect sizes cannot extend to infinity. The power law distribution increases the computational demand of microstructural evolution simulations, due to the necessity to include defects of ever increasing size, making it in practice necessary to truncate the distribution [5], which to date has been done in an ad hoc manner. On the other hand, ion irradiation experiments [3] reveal a clear change in the slope of the distribution in the limit of large defect sizes. The reason for the abrupt decline in the frequency of occurrence of defects above a certain critical size has remained an open question.
In this Letter we develop a model that explains why the statistics of defect cluster sizes deviates from the simple power law behaviour given by Equation (1). The model is fundamentally based on the notion of formation of subcascades, since the cascade-induced formation of large defect clusters in metals occurs in the dense heat spike regions, due to the dynamics of the recrystallization front in the heat spike [6][7][8][9][10]. Hence the formation of such defects depends sensitively on subcascade splitting. It is well known that low-energy (with energies from 1 to 10 keV) recoils in heavy metals produce damage in nearly spherical heat spikes [6,8]. With increasing recoil energy, the cascades start splitting up into spatially separated subcascades [11,12]. This is not a sharp transition: even above the 'subcascade splitting energy' where subcascade formation can be observed, some cascades remain compact, i.e. they still retain the form of a single approximately spherical collision region [13]. Conversely, below the threshold energy, cascades may also occasionally exhibit partial splitting. This gradual transition occurs because the nuclear collision cross section gradually decreases as a function of energy, leading to recoil atoms on average travelling increasing distances between successive collisions [14]. There is nevertheless a non-zero probability for a close succession of collisions to occur, concentrating the cascade energy to a single volume, even when the recoil energy is above the level sufficient for cascade splitting.  Figure 1 illustrates the development of a heat spike and the subsequent defect formation. At the peak of a heat spike, the disordered cascade core is surrounded by an expanding compressed region that deforms the material far outside the liquid core. The denser region shows as a layer of multiply occupied Wigner-Seitz cells surrounding the cascade, with no accompanying vacancies. The particularly large defect clusters remaining after 15 ps nevertheless form along the perimeter of the disordered core region, while the compressed region far outside the core relaxes back to perfect crystal. Hence the disordered, or liquid-like, volume of subcascades defines a natural upper bound on the size of the defects that may form athermally as a result of rapid quenching of the initially overheated molten cascade zone.
Size-frequency distributions of defects derived from MD simulations of cascades in W for various primary knock-on atom (PKA) energies are shown in  Single defects, and clusters containing only a few (N = 2−5) defects, also can be fit to power laws, but with a slope steeper than that exhibited by larger defects, for both vacancies and self-interstitial atoms (SIAs). This is due to different mechanisms governing the formation of very small clusters and point defects [15].
The upper limit of the distribution is determined by the volume of the cascade, or by the volume of the largest subcascade pocket in the event of subcascade splitting. To find the largest possible defect size that can form in a given (sub)cascade, we use the spherical liquid volume approximation. Noting the fact that interstitial clusters in W form two-dimensional dislocation loops, and assuming the loops have the 1/2 111 Burgers vector, with the largest possible loop diameter D equal to the diameter of the liquid volume, we arrive at the maximum possible size in terms of point defects N int D of an interstitial-type defect as [3] where a 0 = 3.165 Å is the lattice parameter of W and D is the diameter of the liquid region. Formation of vacancy clusters occurs inherently differently from SIA clusters, for example the scaling law exponents in Figure 2 are different for vacancies and SIAs. Almost all of the vacancies are formed within the liquid core of a cascade, and are pushed towards the centre by the recrystallization front [9]. In most cases, a single (often diffuse) vacancy cluster is formed in the centre of a spherical cascade volume, usually surrounded by a cloud of single vacancies. Multiple vacancy clusters form only if (partial or complete) subcascade splitting has occurred. In the events where clustering is especially effective, almost all the vacancies condense into one large vacancy cluster. By the law of mass conservation, the largest conceivable vacancy cluster, containing all the vacancies formed in the cascade, must thus contain the same number of vacancies as the total number of individual SIAs produced in the cascade. The SIA clusters form around the surface of the liquid region, and from visual inspection of numerous cascades, we expect that, at most, the SIA clusters may cover half of the surface area, which is equivalent to twice the projected area of the sphere. In addition, we disregard the relatively small numbers of isolated crowdions formed far from the heat spike as a result of replacement collision sequences. Hence the number of vacancies can reach twice the number of point defects contained in the largest possible SIA cluster, and we can relate the diameter D of the liquid region to the largest possible vacancy cluster N vac D according to We point out that this formula does not assume any specific geometry for the vacancy cluster. It is thus valid irrespective of whether the final cluster is a threedimensional 'depleted zone' observed in detailed FIM experiments [16], or a dislocation loop as observed in TEM experiments [3].
To find the frequency f (N D ) of occurrence of a cluster of size N D , we now need to find the probability that the cascade will generate a liquid volume large enough to form such a defect, i.e. a volume larger than or equal to size D. The fractal nature of cascades [17], and the gradual transition from compact cascade volume to subcascade splitting, leads us naturally to look for a power law distribution of subcascade sizes, which in the upper size limit reaches the critical point N c where all the cascade energy is contained in a single compact cascade. This suggests introducing a critical exponent κ combined with a 'reduced size' of the subcascade n = (N c − N)/N c . We thus expect the frequency of formation of a subcascade large enough to accommodate a defect of size N to be given by where F SD (N) is the frequency of occurrence of a subcascade of size D(N) (where D is related to N through Equation (2)), and B is the average total number of subcascades. In the limit of small N, every subcascade is large enough, and accordingly lim N→0 f SC = B. The frequency of occurrence of a defect of size N, in cascades initiated by a given PKA energy above the subcascade splitting threshold, can now be weighed by the probability that a subcascade of sufficient size will form, giving the total size-frequency distribution of defects as Here we choose not to make any generalizing assumptions, but rather consider the parameters as material-and PKA energy-dependent, and determine them by fitting to simulation data. We find parameters for the distribution of subcascade sizes (Equation (4)) using a recently developed method, detailed elsewhere [18], based on the binary collision approximation (BCA) [19]. We have simulated 200 cascades with the BCA method for each PKA energy. The frequency of occurrence of subcascades larger than N predicted by these calculations is shown in Figure 3, and follows very closely the distribution in Equation (4).
As the PKA energy decreases below the subcascade splitting threshold, the distribution approaches that of a Heaviside step function rather than following Equation (4), as can be seen in Figure 3 for 100 keV cascades. The BCA subcascade distributions are the same for SIA and vacancy clusters, with only the N c parameter changing due to the different mechanisms of cluster formation, leading to different maximum cluster sizes. The final distributions are different for vacancies and SIAs due to the different scaling laws derived from MD simulations ( Figure 2). To find the final distributions, we assume that the MD results are accurate in the defect size range that they cover. Furthermore, we recognize that partial subcascade splitting already affects the distributions derived from MD at these PKA energies. We therefore fit the parameters S and A in the full distribution function (Equation (5)) to the MD results, with parameters B, N c and κ previously determined by the subcascade distribution. We recall that the frequency of defects in MD simulations is proportional to the PKA energy, with the same slope (Figure 2), and thus multiply the frequency in 200 keV cascades by 2, to obtain the distribution for the 400 keV cascades. The final parameters for the best fit are given in Table 1.
In bulk cascades, the total numbers of vacancies and SIAs should be equal, despite the different statistics. As a consistency check, we calculate the total number of defects N tot of each type, by integrating Nf (N)dN over the full range of the distributions, where f (N) is given by Equation (5). Since point defects and small clusters in the range N ∼ 2-3 follow different scaling laws, we integrate over 4 ≤ N ≤ N c , and use the clustered fraction of defects F Cl derived from MD simulations to determine Notes: The parameters B, N c and κ describe the subcascade splitting, and were determined from the BCA data alone, after which A and S, describing the defect sizes within (sub)cascade regions, were determined by fitting the complete model in Equation (5) to the MD data. The error bars are statistical uncertainties obtained using the Levenberg-Marquardt least-squares fitting algorithm [20]. Fitting was performed assuming a 10 % statistical uncertainty for all the data points. Figure 4. The size distribution of defects as predicted by MD simulations of irradiated thin foil, combined with BCA simulations in the full model presented in this paper. For comparison, also the predictions of a simple power law, fitted to the same MD data, is given. Lines labelled 'with loop loss correction' show the distributions fitted to the same MD and BCA data, with the added correction for loss of mobile loops to the surface which occurs on longer time scales than those of the MD simulations, during the time taken for TEM experiments (see text for details). The experimental results of self-ion irradiated W foils are also shown, but has not been used in the fitting of the lines. Error bars for the experimental data express the uncertainty in the loop count from the micrographs (details given in the text). The size in nm is indicative. The vertical line marks the experimental limit of visibility [21].
N tot . The results are given in Table 1. Good agreement is found between the predictions for SIAs and vacancies, although a relatively large uncertainty in the fraction of clustered vacancies in 150 keV cascades gives rise to a large uncertainty in the total vacancy count. Figure 4 shows the predictions of our model for dislocation loops produced in 150 keV cascades, compared to TEM observations of self-ion irradiated W. Irradiations were performed in situ on the IVEM-Tandem Facility at Argonne National Laboratory, with 150 keV W + ions at 30 K, up to 1.25 × 10 16 W + /m 2 (equivalent to about 0.01 dpa [22]) at around 15 degrees off the zone axis of ultra-high-purity, well-annealed (001) grains. Details of the experimental technique can be found in [3]. Each foil region was imaged using weak-beam dark-field conditions, using (g = 200, 4.25g; g = 200, 4.75g; g = 200, 5.25g; g = 110, 5.25g; g = 110, 7.25g; and g = 110, 7.75g).
The regions were superimposed using the convergent weak beam technique of Prokhodtseva et al. [23] prior to automated loop counting using the method described in [3]. Eight regions of the foil were studied, with a combined area of 3 µm 2 , and a total of 37,500 incident ions in the observed regions. The error bars of the experimental observations in Figure 4 are derived from a procedure which corrects for miscounting due to the faintness of spots near the visibility limit (see Supplementary).
To account for the surface effects, the parameters A and S of Equation (5) were refitted to MD simulations of foil irradiation reported in [15], performed with the same geometry as in the experimental set-up, using a simulation cell with an open surface on top and bottom, with grain and ion orientation as in the experiments, and a depth of the simulation cell comparable to the thickness of the TEM foils. The resulting distribution given by Equation (5), with A = 4.01 ± 0.62 and S = 1.78 ± 0.04, is plotted in Figure 4, together with the MD data points from [15], and the simple power law distribution that fits the MD data.
The MD simulations of foil irradiation account for the effect of the surface on the formation of primary defects, but the time scale of the simulations is too short to account for the loss of mobile loops to the surface which occurs during the time it takes to produce the TEM micrographs. We account for the longer time scale of experiments by correcting the prediction for the loss of mobile loops to nearby surfaces. Based on the spatial distribution of defects in the MD foil simulations, we assume an initial depth of the loops given by a Gaussian probability distribution centred at d = 15 nm depth, with a half-width of 10 nm. The average diameter of a 150 keV cascade is 10 nm. We determine the number of loops remaining in the foil based on the probability that a loop will be elastically trapped given the degree of thermal fluctuations at 30 K, and an observation time of t = 10 min. The corrected prediction is also plotted in Figure 4, with the full calculation given in the supplementary information. For simulations of microstructural evolution under irradiation, the uncorrected distribution naturally should be used, since the correction in effect accounts for the evolution, and in fact gives the same results as OKMC simulations building on the same assumptions of loop mobility and trapping [24].
From the form of the function describing the probability of a loop being elastically trapped in the foil (Equation (8) of the supplementary information), it can be seen that the self-trapping probability is most strongly dependent on the size and count of loops, followed by the elastic constants and temperature, but is only weakly dependent on the observation time. We have also found that the self-trapping probability is only weakly dependent on the size of the subcascade (σ ), cascade depth, or on whether the defect distribution is taken to be Gaussian (as is done here) or homogeneous. Hence the corrected distribution is most strongly dependent on the parameters determined from the MD data presented here, while variations, within reasonable bounds, of the other parameters used in the correction function have negligible effect on the final curve. No fitting to experimental data was required to achieve the good match between observed and predicted loop count in Figure 4.
The experimental observation limit, determined by the number of TEM foils and the total fluence, is also indicated in Figure 4 as the limit of one loop seen per bin, showing that the experimental data is able to distinguish between the predictions of a power law and the full model, and agreement is found with the full model. The simple power law overpredicts defect numbers in the range of large defects, which is a critical part of the size distribution, since the larger clusters are the ones most likely to be trapped in the foil, thus contributing the most to the accumulating damage.
In conclusion, we have proposed and validated a model for the distributions of sizes of SIA and vacancytype defects produced by energetic cascades in irradiated metals. We have parameterized the model based on MD and BCA simulations for the specific case of tungsten. The model successfully extends the predictions of defect sizes beyond the N < 200 range directly accessible by the MD simulations. In addition, the versatility of the model, derived for bulk simulations, is demonstrated by direct comparison to foil irradiation experiments, where the effects of the surface on the primary damage formation are accounted for by fitting to relevant MD simulations. By further introducing a probabilistic correction for the loss of mobile loops to the surface over the longer time scales of experiment, we achieve very good agreement with TEM observations of self-ion irradiated tungsten foils. The model is applicable to high PKA energies above the subcascade splitting threshold, in particular for energies not easily accessible to direct MD simulations. The distributions are different for different types of defects, owing to the difference between the dynamics of their formation. The model provides well-defined upper limits on defect sizes occurring in cascades, which agrees with experiment and is critically important for the multiscale modelling of microstructural evolution of materials under neutron and ion irradiation.