In-cascade recombination effects on defect production

The physical basis for the overprediction by the NRT-dpa model of the defect production at high energies is the enhanced recombination of defects in close proximity in energetic displacement cascades. The binary collision simulations used as the basis of the NRT-dpa model3 focused on the collisional phase of the displacement cascade and did not consider the dynamics of cascade evolution as atomic velocities fell to the speed of sound (~5 eV) and lower, when many-body interactions become relevant. In energetically dense cascades, like that shown in Fig. 1, local melting clearly plays an important role in defect retention and structure. With increasing primary knock-on atom energy, the displacement event produces progressively more Frenkel defects (pairs of vacancies and interstitials25) that are spatially close to other defects. The ~10–100 jumps occurring per atom during the 1–10 ps cascade cooling phase14 can induce significant additional recombination events as the cascade atom energies decrease, following the collisional phase, from E d to the threshold value for atomic migration (E m ~ 0.01–0.3 eV for self-interstitial atoms and ~0.5–1 eV for vacancies). Accurate simulations of these cooperative multi-body effects in displacement cascades are realistically performed with MD simulations14.

Figure 2a summarizes the defect production as a function of primary knock-on atom energy as determined from experiments performed in Cu near 4 K (where long-range thermally activated defect motion is impossible25). The predicted defect production and number of replaced atoms obtained from MD simulations are also shown. The figure shows that the actual defect production is sublinear with respect to damage energy between ~0.1 and 10 keV (ref.23), becoming about 1/3 of the NRT-dpa prediction. At energies >10 keV corresponding to the onset of subcascade formation14,26,27, the defect production increases linearly with damage energy but maintains the factor of ~3 lower defect production compared to the NRT-dpa value.

The physical basis for the reduction in surviving defects, with respect to the NRT model, with increasing knock-on atom energy can be understood by considering the following simplified derivation.

The ultimate survival of initially created Frenkel defects requires physical separation of the interstitial and vacancy beyond a minimum distance known as the spontaneous recombination distance (L). Atomic collisions along close-packed directions (known as recoil collision sequences) are one example of a method to efficiently transport interstitial atoms to the periphery of a displacement cascade, leaving the associated vacancy near the cascade interior. Molecular dynamics simulations22 indicate atom transport from the displacement cascade interior may be associated with a supersonic shock-front expanding from the primary recoil event during the early stages of the cascade evolution. At low energies (below the subcascade formation regime14) the displacement cascades are roughly spherical with radius R, and form a liquid-like zone of dense collisions (the heat spike described above).

It is further assumed that only interstitials transported to the cascade outer periphery defined by R − L to R will result in stable defects, whereas Frenkel pairs created in the cascade interior (0 to R − L) will experience recombination. The fraction of initially created NRT-dpa defects that survive is therefore given by the ratio of the outer spherical shell volume to the total cascade volume:

$$\begin{array}{l}\xi _{{\mathrm{survive}}} = \frac{{V_{{\mathrm{outer}}} - V_{{\mathrm{inner}}}}}{{V_{{\mathrm{outer}}}}} = \frac{{{{(4{\mathrm{\pi }}R^3}}/{3}) - {{((4{\mathrm{\pi }}(R - L)^3}})/{3})}}{{{{4{\mathrm{\pi }}R^3}}/{3}}}\\ = 3\frac{L}{R} - 3\frac{{L^2}}{{R^2}} + \frac{{L^3}}{{R^3}} \approx 3\frac{L}{R}\end{array}$$ (2)

for L ≪ R. This ‘surviving defect production fraction’ \(\xi _{{\mathrm{survive}}}\) thus tells which fraction of defects predicted by the NRT-dpa model without any recombination survives. The cascade radius R can be, within the regime of spherical cascades, estimated from classical theory of nuclear stopping power28,29. In practice, we used the SRIM code that implements an integral calculation to obtain mean range tables, based on cross sections from the widely used Ziegler−Biersack−Littmark (ZBL) interatomic potential29.

We found that low-energy (less than or of the order to 10 keV) recoils of damage energy T d have an average movement distance (range) R that is proportional to \(T_{\rm d}^x\), where the exponent x is ~ 0.4–0.6 for the metals considered in this study. Since \(R\propto T_{\mathrm{d}}^x\), this further gives

$$N_{\rm d} \prime \left( {T_{\rm d}} \right)\frac{{0.8T_{\rm d}}}{{2E_{\rm d}}}\xi _{{\mathrm{survive}}} = \frac{{0.8T_{\rm d}}}{{2E_{\rm d}}}3\frac{L}{R}\propto \frac{{0.8T_{\rm d}}}{{2E_{\rm d}}}3\frac{L}{{T_{\mathrm{d}}^x}} \propto T_{\mathrm{d}}^{1 - x}.$$ (3)

This simple model thus provides an intuitive explanation for why cascade damage production is sublinear with damage energy in the heat spike regime. Physically realistic MD simulation studies14,30 have reported that defect production rates up to the onset of subcascade formation in a variety of metals can be well described by N d ~ (T d )1−x, where x is between 0.2 and 0.3. These x values are slightly larger than the value obtained in our simplified model because real cascades are not perfectly spherical and some defects form small clusters, reducing the recombination probability.

However, it is well known that at high energies cascades break up into subcascades24,31,32, after which damage production becomes linear with energy. Hence the surviving defect fraction factor \(\xi \left( {T_{\mathrm{d}}} \right)\) that accounts also for subcascade breakdown should have the feature of being a power law at low energies, but becoming a constant c at high ones. A function that fulfils both criteria is

$$\xi \left( {T_{\mathrm{d}}} \right) = A\prime T_{\mathrm{d}}^b + c,$$ (4)

where b < 0 is consistent with the damage production efficiency reducing with increasing energy T d and the desired limit \(\xi \left( {T_{\mathrm{d}}} \right) \to c\) when T d → ∞. This thus gives a total damage production

$$N_{\mathrm{d}}\prime \left( {T_{\mathrm{d}}} \right) = \frac{{0.8T_{\mathrm{d}}}}{{2E_{\mathrm{d}}}}\left( {A\prime \,T_{\mathrm{d}}^b + c} \right) = \frac{{0.8A\prime \,T_{\mathrm{d}}^{1 + b}}}{{2E_{\mathrm{d}}}} + \frac{{0.8\,c\,T_{\mathrm{d}}}}{{2E_{\mathrm{d}}}}.$$ (5)

Note that here the exponent b is not the same as x, since the latter ξ function is not a pure power law. The prefactor A′ is defined by demanding the function to be continuous, i.e. \(\xi (2E_{\mathrm{d}}/0.8) = 1\).

Taken together, this derivation leads us to propose (based in part on review work done within an OECD Nuclear Energy Agency group24) a modified defect production model, the athermal recombination corrected displacements per atom (arc-dpa).

$$N_{{\mathrm{d}},{\mathrm{arcdpa}}}\left( {T_{\mathrm{d}}} \right) = \left[ {\begin{array}{*{20}{c}} 0 & , & {T_{\mathrm{d}} < E_{\mathrm{d}}} \\ 1 & , & {E_{\mathrm{d}} < T_{\mathrm{d}} < \frac{{2E_{\mathrm{d}}}}{{0.8}}} \\ {\frac{{0.8T_{\mathrm{d}}}}{{2E_{\mathrm{d}}}}\xi _{{\mathrm{arcdpa}}}\left( {T_{\mathrm{d}}} \right)} & , & {\frac{{2E_{\mathrm{d}}}}{{0.8}} < T_{\mathrm{d}} < \infty } \end{array}} \right]$$ (6)

with the new efficiency function \(\xi _{{\mathrm{arcdpa}}}\left( {T_{\mathrm{d}}} \right)\) given by

$$\xi _{{\mathrm{arcdpa}}}\left( {T_{\mathrm{d}}} \right) = \frac{{1 - c_{{\mathrm{arcdpa}}}}}{{\left( {2E_{\mathrm{d}}{\mathrm{/}}0.8} \right)^{b_{{\mathrm{arcdpa}}}}}}\,T_{\mathrm{d}}^{b_{{\mathrm{arcdpa}}}} + c_{{\mathrm{arcdpa}}}.$$ (7)

Here E d is the average threshold displacement energy33 which is the same as in the NRT-dpa and b arcdpa and c arcdpa are material constants, that need to be determined for a given material from MD simulations or experiments. The overall form (Eq. (5)) and the constant 0.8 are retained for direct comparison with the NRT-dpa model; in particular making it easy to modify computer codes that now use the NRT-dpa by simply multiplying with the function \(\xi _{{\mathrm{arcdpa}}}\left( {T_{\mathrm{d}}} \right)\).

Figure 3 compares the derived arc-dpa expression for Fe and W with several recent MD simulation results used for the fitting. We tested that if the fit is limited to energies <10 keV, one also can fit the data well with a power law with an exponent of ~0.7–0.8, i.e. the data is consistent with MD reports of power law dependencies. However, the arc-dpa form has the major improvement that it can also describe the saturation. Even though there is some variation in the MD data (due to differences in interatomic potentials), all of the MD results give damage production well below the ξ = 1 value predicted by the NRT-dpa model for cascade energies >1 keV. The arc-dpa fit to the composite data gives a reasonable averaging description of the decreasing trend in ξ up to ~10 keV and the expected approach to a constant value at higher cascade energies.

Fig. 3 Improvement with arc-dpa and rpa. Illustration of the improvement obtained with the new arc-dpa and rpa equations for a Fe and b W. The W data also includes two data points simulated at 800 K with the DD potential (solid circles). The references are: [A98]: ref. 26, [Z93]: ref. 13. The Fe damage data is from ref. 14 (Stoller) and ref. 48 (AMS, MEA-BN, DD-BN). The Fe replacement data and all W data is original work for this publication, see Methods section Full size image

Replacements-per-atom (rpa) model

Since the NRT-dpa model deals with production of defects that are not on perfect lattice sites, it cannot predict the number of atoms that are transported from their initial lattice site to a new lattice site, i.e. replace another atom in a perfect crystal site (right panel in Fig. 2c). This number of atom replacements is experimentally measurable via so-called radiation mixing experiments. Typically, an ion beam is used to bombard a thin marker layer inside a material, and the resulting broadening of the marker layer is measured16,34. For ordered alloys, it can also be conveniently measured by electrical resistivity15. Via an analogy with random walk atom diffusion, it is possible to relate this measured broadening to the actual number of atom replacements per ion inside the material35. Analysis of neutron and ion beam radiation mixing data has shown that the actual number of replaced atoms can be more than an order of magnitude larger than the number of displacements predicted by the NRT-dpa model15,36,37,38. A correct estimation of this number can be of enormous importance in predicting the effects of irradiation on phase stability39,40 and the associated mechanical properties of materials. Nanostructured materials, such as nanolaminates and nanoscale oxide-dispersion strengthened steels, are particularly sensitive to these errors owing to their small length scales.

The superlinear increase in the number of replaced atoms with increasing knock-on atom energy can be understood by a model considering the spatial extent of a collision cascade. We consider first low energies (in the keV regime) and dense materials, where cascades are normally compact. As noted above, low-energy cascades are roughly spherical. After the ballistic phase of a cascade, MD simulations show (cf. Fig. 1) that the lattice breaks down and a liquid-like region forms. In this region all atoms are free to move and hence are almost certain to lead to one or more replacements during the thermal spike phase (as illustrated in Fig. 2c, right frame). The number of atoms N in a spherical cascade of radius R is proportional to the sphere volume, i.e. N∝R3, and (as already noted for the arc-dpa model) \(R \propto T_{\mathrm{d}}^x\). We thus find that the number of replaced atoms \(N_{{\mathrm{rpa}}} \propto T_{\mathrm{d}}^{3x}\). Since x > 1/3, this simple consideration gives an intuitive explanation for why the number of replaced atoms increases superlinearly with energy at low energies, when cascades are compact. At high energies, when cascades split into subcascades32,41, the behaviour can be expected to change to a linear dependence with damage energy. Similarly to the arc-dpa function, we thus arrive at a form

$$\xi _{{\mathrm{rpa}}}\left( {T_{\mathrm{d}}} \right) \propto \frac{{T_{\mathrm{d}}^{c_{{\mathrm{rpa}}}}}}{{b_{{\mathrm{rpa}}}^{c_{{\mathrm{rpa}}}} + T_{\mathrm{d}}^{c_{{\mathrm{rpa}}}}}}.$$ (8)

The proportionality prefactor is again set to ensure continuity, (2E d /0.8) = 1. To augment the NRT equation in order to predict the number of replaced (mixed) atoms, we propose another improved damage function, the replacements-per-atom (rpa) equation. The correction factor is applied in Eq. (5) as for the arc-dpa, but now it has the form

$$\xi _{{\rm rpa}}\left( {T_{\rm d}} \right) = \left( {\frac{{b_{{\rm rpa}}^{c_{{\rm rpa}}}}}{{\left( {2E_{\rm d}{\mathrm{/0}}{\mathrm{.8}}} \right)^{c_{{\mathrm{rpa}}}}}} + 1} \right)\frac{{T_{\mathrm{d}}^{c_{{\mathrm{rpa}}}}}}{{b_{{\mathrm{rpa}}}^{c_{{\mathrm{rpa}}}} + T_{\mathrm{d}}^{c_{{\mathrm{rpa}}}}}}.$$ (9)