Cluster dynamics simulation of deuterium retention behaviors in irradiated beryllium

Abstract

The retention behaviors of deuterium (D) in beryllium (Be) are investigated using a spatially resolved cluster dynamics model under different irradiation conditions. The trapping effects of deuterium (D) in the forms of D atoms, D₂ molecules and D with vacancy complex clusters (D_mV) play the most important role in the behaviors of D retention in bulk Be under irradiation of 9 keV D ions. The fraction of D₂ in the total D retention increases with the increase in ion influence, due to the chemical reaction rate enhancement between D atoms with high density. The increases in both ion incident angle and Be bulk temperature reduce the retention of D_mV complex clusters by increasing the D desorption rate. In addition, the neutron synergistic irradiation changes the D retention profiles, especially in the recombination region, by introducing extra defect sinks. These results can improve the understanding of the mechanisms of D diffusion, accumulation and retention in irradiated Be.

---

1 Introduction

Nowadays, a shortage of energy and environmental issues have prompted people to spend great effort on the research of nuclear fusion energy. Energetic particles, such as hydrogen isotopes and helium ions, that escape from the fuel during plasma burst events in fusion reactors, would implant into and damage the plasma-facing materials (PFMs).^1,2 The aggregation of different radiation defects, especially incident ion deposition, can induce strong disruption in the host lattice, thus leading to degradation of PFMs.³ Beryllium (Be) is considered as one of the promising candidates for PFMs in nuclear fusion reactors, because of its excellent oxygen gettering capability, low sputtering, great mechanical and thermal properties, etc.⁴ However the high retention rate of deuterium (D) in Be induces a set of serious problems, such as the formation of surface sputtering, erosion, swelling, blistering, etc. To understand the D retention behaviors in Be, it is necessary to predict the detailed interactions of D with other defect clusters in Be in a quantitative manner.

Numerous experimental studies have been carried out to reveal the accumulation behaviors of radiation damage. For instance, the Elastic Recoil Detection (ERD) method has been used to study the D retention behaviors in bulk Be,⁵ which shows that D can penetrate to a depth of 450 nm and accumulate up to a concentration of (7–9) × 10²¹ D cm⁻³ under a D flux of 10¹⁶ D cm⁻² s⁻¹. After that, similar experiments illustrated that deuterium is mainly trapped as D atoms and D₂ molecules in irradiated Be.⁶ On the other hand, theoretical studies have also been conducted, focusing on the atomic details of D behaviors, such as the solubility, migration paths and correlative kinetic coefficients of D in Be, by using atomic methods such as density functional theory (DFT) and molecular dynamics (MD).^7–11 However, the long-term behaviors of D in Be and the details of dynamical evolution are inadequately studied, which requires the employment of a multi-scale model.

In materials science, the cluster dynamics (CD) model is a mesoscopic model to describe the evolution of defect concentrations by considering the generation, diffusion, reaction and absorption processes of point defects and clusters with a possible event list and corresponding rate coefficients. Compared to the atomic-scale simulation methods, DFT and MD methods, the CD model is more effective for handling the long-term irradiation damage evolution, which encompasses models from the atomistic to the continuum scale and timescales from diffusion (∼ns) to micro-structural evolution (∼years). The high computational efficiency is due to the introduction of the mean field approximation, which refers to a basic hypothesis of uniform distributions of defects (the real system is replaced by an effectively continuous medium), but with the loss of the spatial correlation effect between defects. We can also partly take into account the spatial correlation effect by partitioning the depth into small intervals in which the defect concentration changes through diffusion and reactions. Additionally, the master equations are used to describe the concentration evolution of different point defects and their clusters with time and space. By numerically solving the set of partial differential equations (PDEs) constituted by the master equations, information about the different defect evolutions can be obtained and compared with the experimental results straightforwardly. The method has been successfully applied to simulate the long-term evolution of the microstructures of materials.^12–17 Thus, it is very suitable to use this method to study the defect dynamics evolution in irradiated systems.

In this paper, we intend to investigate the long-term evolution of D diffusion and retention with depth for Be under different irradiation conditions, by employing the CD model. The model is extended from the IRadMat program^12,16,18 by taking into account the corresponding types of defects and by adopting reliable parameters and complex reaction types. Our results highlight the different roles of ion fluence, incident angle and system temperature as well as the neutron synergistic effect on the D retention behaviors in Be.

2 Model and method

We adopt a deterministic cluster dynamics model based on the mean-field rate theory to simulate the retention behaviors of D in Be under energetic particle irradiation.^12,13 A set of partial differential equations (PDEs) constituted by the master equations of different defects in different volume units must be solved numerically. The defect distributions with depth obtained using the model can be compared with experimental results, which could provide information about the roles of the different reaction mechanisms. To simplify the calculation, we assumed that the self-interstitial atom (SIA, I), di-interstitial (I₂), vacancy (V), deuterium (D) and di-deuterium (D₂) defects are mobile, while the complex defect clusters of I_n (n ≤ 100), DI and D_mV (m ≤ 5) are immobile. Additionally, we only consider the defect distribution with depth for two main reasons. On the one hand, handling a 3D irradiation system with a deterministic CD model will increase the computational cost significantly. On the other hand, since the Be surface is irradiated under a uniform distribution D ion beam in the experiment, the total D fluence and the induced defects will not change along with the planes parallel to the surface. Thus, it is a very reasonable approximation to reduce a 3D model to a 1D model for simplification, as commonly employed elsewhere.^13,16 It is reasonable to neglect the minor intrinsic thermal defects at the beginning of the simulation and consider the first-order boundary (corresponding to C (z = 0) = 0) for simplicity because of the relatively high diffusivity of D in Be. In addition, the diffusion coefficients are anisotropic along the different specific directions in a hexagonal structural system like Be and Zr, which have been studied by Woo et al. and Barbu et al.^19–21 However, it is noteworthy that this anisotropy is only considerable for specific defects, especially SIAs, with extremely high diffusivity.²⁰ In addition, the contribution of SIAs to D retention in Be can be neglected, as shown below. Thus, we have ignored the influence of the defect anisotropic diffusivity for the D retention in Be. The basic procedure is demonstrated briefly as follows.

2.1 Master equation

The master equation describes the evolution of defect concentrations with time in the irradiated system, including the generation, diffusion, reaction and absorption of different kinds of defects,^{12,13,16,17,22} which is given by,where C_θ (θ = I, I₂, V, D, D₂, I_n, DI, D_mV) is the concentration of defects θ in the irradiated system at a specific depth and time; G_θ, D_θ, and L_θ represent the production rate, diffusion and inherent absorption coefficients (including dislocations and grain boundary) of the defects θ, respectively. w(θ, θ′) is the transition rate coefficient per unit concentration of defect cluster θ into defect cluster θ′. In order to further decrease the computational cost, the Fokker–Plank approximation is adopted in our model to transform these discrete master equations into continuous equations.

2.2 Rate coefficients

According to the above assumptions, the reaction types and rate coefficients are listed in Table 1. The rate coefficients describe the probability of occurrence of the corresponding reactions. The forward and backward chemical reaction rate-coefficients can be calculated using the theory of diffusion-limited reaction and chemical equilibrium principles, respectively. In addition, the generation rate of point defects during irradiation, G_I/V, is determined by using TRIM code.²⁸

Table 1 Reaction types and the relevant chemical rate coefficients. Reaction coefficients quantify the rate of a chemical reaction for unit concentration, which are mainly dependent on system temperature (T), atomic volume (Ω), cluster binding energy (E_b), defect relative diffusivity (D) and reaction distance (R_A+B). In principal, the reactions are described by A + B ⇌ AB with the forward and backward rate coefficients calculated from k_A+B⁺ = 4πR_A+B(D_A + D_B) and , respectively

Reaction types	Coefficients rate
I + V ⇌ 0	k^+I+V, GI/V
I + In ⇌ In+1	α^+n, αn+1^−
I + D ⇌ DI	k^+D+I, kDI^−
V + D ⇌ DV	k^+D+V, kDV^−
I2 + V → I
D + D ⇌ D2	η^+, η^−
I + DmV → mD
I2 + In ⇌ In+2	β^+n, β^−n+2
V + In ⇌ In−1
V + DI → D	k^+DI+V
D + DmV ⇌ Dm+1V	γ^+m, γ^−m+1
I + L → LI	LI
I2 + L → LI2	LI2
D + L → LD	LD
D2 + L → LD2	LD2

---

2.3 Numerical method

The master equations describing the evolution of different defects compose a set of PDEs. To improve the computational efficiency, the PDEs can be transformed into a set of ordinary differential equations (ODEs) by the Taylor series expansion up to second-order terms.^29,30 In this paper, the set of ODEs is solved by using lsoda subroutine packages.³¹ The method is efficient enough for handling ∼10³ ODEs here, with no more than one hour on a modern personal computer.

3 Results and discussion

The CD model is a mesoscopic method,³² in which most of the input parameters should be suitably chosen from the values of experiments or atomic calculations, and the recommended parameters are listed in Table 2. The accuracy of the CD model is mainly dependent on the reliability of the parameters selected, especially the critical parameters including the defect formation energies and migration energies of mobile point defects (I, I₂, V, D and D₂) which dominate the reaction dynamics processes. By considering the published atomic-scale computational and experimental results, the critical defect types, the reaction event lists and the corresponding reaction coefficients are carefully selected. For example, in order to fit the experiments very well, the D migration energies (E^m_D) are adjusted in the reasonable range of about 0.29–0.41 eV and set as 0.41 eV,^4,8,27 by considering that the presence of impurities and inherent defects probably affects the incident D diffusion in Be, as pointed out by Ortiz et al.³³ In addition, these parameters have been commonly adopted elsewhere.³⁴ The ab initio results predicted that the maximum number of D surrounding a vacancy to form a stable cluster is five with binding energies between 1.27 and 0.42 eV, as listed in Table 2.⁸ In addition, the binding energies of mobile point defects (I, I₂) with different types of large loops (I_n) can be obtained using the capillary law approximation.^35,36

Table 2 Correlated parameters used in the case of D ions and neutrons irradiated on Be

Parameters	Symbol	Value	Ref.
D beam intensity	ID	10^18 to 10^20 m^−2 s^−1	21
Temperature	T	300 K	23
Lattice parameter	a0	2.27 Å	24
	c0	3.56 Å
D radius	rD	0.53 Å	13
Burgers vector	b	1.78 Å	24
Dislocation line density	ρD	10^13 m^−2	25
Recombination radius	rIV	4.54 Å	—
Formation energy of SIA	E^fI	5.24 eV	26
Formation energy of vacancy	E^fV	0.81 eV	27
Formation energy of D	E^fD	1.71 eV	27
Migration energy of SIA	E^mI	0.02 eV	27
Migration energy of vacancy	E^mV	0.8 eV	27
Migration energy of D	E^mD	0.41 eV	4
Binding energy of D2		2.3 eV	27
Binding energy of D–V	E^bDV	1.27 eV	8
Binding energy of D–D2V		0.93 eV	8
Binding energy of D–D3V		0.77 eV	8
Binding energy of D–D4V		0.54 eV	8
Binding energy of D–D5V		0.42 eV	8

---

To verify our model, we plot the D retention profile with implant fluence and the concentration profile of D trapped as D atoms and D₂ molecules in the near-surface region to compare with experiments. It is shown in Fig. 1(a) that the total amount of deuterium retention increases linearly with increasing D fluence, which is in good agreement with the experimental results when the D beam fluence is low.^23,37 However, for high fluence, the experimental relation deviates from linear to reach a constant value, due to the surface reconstruction under high flux D ion erosion which prevents the diffusion of implanted D into bulk regions. Furthermore, the saturation of the trapping sites also makes the retention decrease when a much higher density of D atoms is accumulated in the near-surface region.^13,37

Fig. 1 (a) The total amount of deuterium and fraction quantities of deuterium trapped as D atoms and D₂ molecules. (b) Concentration profiles of deuterium trapped as D atoms and D₂ molecules in a flux of 3 × 10¹⁹ D m⁻² s⁻¹ up to a fluence of 1.9 × 10²¹ D m⁻². D atoms are implanted into beryllium with an energy of 9 keV at 300 K.

Fig. 1(b) shows the concentration of deuterium against depth near-surface. Deuterium prefers to be in the form of D₂ molecules trapped in the near-surface of Be under D ion irradiation at a flux of 3 × 10¹⁹ m⁻² s⁻¹ up to a fluence of 1.9 × 10²¹ D m⁻². These results compare well with those of experiments,³⁷ except for a discrepancy near the surface region. The enhancement of implanted particles at the near-surface of an irradiation target in measurements is a common phenomenon.¹³ The difference between the experimental and calculated D concentration at the near-surface arises mainly from the surface roughness, which enhances the adsorption of D atoms on the Be surface in the experiments, and partly from the choice of the free boundary condition on the system surface in the calculation. Furthermore, the plasma–surface interaction (PSI) induces surface configuration reconstruction into a new fluctuated surface, which can hamper the implanted D atoms and cause the enhancement of D accumulation at the near-surface. Due to the materials used for nuclear reactors generally being subjected to the generation of a non-equilibrium concentration of atomic defects,³⁸ the concentration ratio between D and D₂ will continually change with different conditions (ion flux, fluence, incident angle, target temperature, etc.) and has a value of roughly 0.5 under the experimental conditions given here.

The distribution of the deuterium with depth in Be irradiated at different D fluences is shown in Fig. 2. With the increase in D fluence, the D deposition extends into a deeper region. There is a peak at several hundred nm, which comes from self-accumulation and trapping by other D-related defect clusters, that is, D_mV, DI and LD. Consequently, the concentrations decrease with increasing depth, which is mainly due to two competitive processes, i.e. the compensation by diffusing from the surface and the loss by diffusing away to the bulk. For a low D fluence of 3 × 10²⁰ m⁻², the D atom concentration is always higher than that of D₂ molecules. However, for high D fluence, the concentrations of D and D₂ reverse at a certain depth, because of the relatively high ratio of D to D₂ concentration, increasing the conversion rate of D + D → D₂.

Fig. 2 The depth profiles of D atoms and D₂ molecules at different fluences under a flux of 3 × 10¹⁹ D m² s⁻¹ with incident energy of 9 keV at 300 K.

To clarify the contributions of the different kinds of D-related defect clusters to the D total depth distribution, we plot the detailed depth profiles of the D-related defect clusters retained in Be in Fig. 3(a), which shows that D atoms, D₂ molecules and D_mV clusters govern the depth distribution of D in Be. The concentration of deuterium trapped in these three forms is over one order higher than those trapped by dislocation lines and SIAs.

Fig. 3 Details of the depth profiles of D retained in Be under a flux of 3 × 10¹⁹ D m⁻² s⁻¹ up to a fluence of 1.9 × 10²¹ D m⁻² with incident energy of 9 keV at 300 K, (a) without and (b) with synergistic neutron irradiation at a damage rate of 10⁻³ dpa per s.

In nuclear fusion reactors, synergistic neutron irradiation will further aggravate damage and hydrogen isotope retention in Be. Its kinetic energy is so high (14.1 MeV) that the neutron in the fusion reactions can easily pass through the PFMs. Consequently, the defects induced by the neutron can be considered as uniform everywhere at a constant rate (G_neutron), which is usually described by the unit of dpa (displacement per atom). We set the neutron irradiation as increasing from a non-zero rate in the non-neutron case to generate defects. Thus, we have also considered this synergistic effect of D and neutron irradiation in the model. The evolution and depth-distribution of defects have been obviously changed under the synergistic neutron irradiation, as shown in Fig. 3(b). The extra point defects induced by neutrons can further trap D atoms, and can compensate for the valley (in the range of several microns) of D_mV clusters formed by the recombination of SIAs and vacancies, as shown in Fig. 3(a). However, the influence of synergistic neutron irradiation on the depth-distribution profiles of the other D-related defect clusters is feeble, due to the weak sink strength. Therefore, the synergetic effect of neutron irradiation can change the D distribution forms considerably, but only slightly the amount of D retention.

Fig. 4(a) illustrates that the distribution of D_mV with depth can be tentatively divided into three regions according to their respective features (see the dashed line divided areas), i.e. (1) a peak in the surface layer (Region I), (2) a valley of several microns (Region II), which can be called the recombination region and (3) a decay with depth entering into the bulk (Region III).^13,39 In Region I, D₅V dominates the concentration of D retained in Be, because of the super-saturation of D in this region during high flux D implantations. In contrast, D_mV clusters prefer to form smaller-size clusters in Region III because of the lower ratio of the concentration of D to V. However, in Range II, almost all of the vacancies are combined by numerous SIAs migrating rapidly from the surface layer, which causes the formation of a sink valley.³⁹ On the contrary, the extra vacancies introduced by synergistic neutron irradiation can compensate for the sink valley, increasing the D retention in this region as shown in Fig. 4(b).

Fig. 4 The depth profiles of deuterium–vacancy complexes D_mV under a flux of 3 × 10¹⁹ D m⁻² s⁻¹ up to a fluence of 1.9 × 10²¹ D m⁻² with incident energy of 9 keV at 300 K. (a) Without and (b) with synergistic neutron irradiation at a damage rate of 10⁻³ dpa per s.

On the other hand, Be usually serves under the conditions of D plasma with random incident angles and different localized temperatures under high-heat loads in nuclear fusion reactors. Since the 1D semi-infinite diffusion-reaction system with free boundary conditions was considered, the D atoms diffuse towards the surface and depart from the bulk directly. In the following, we will consider these two factors for the illustration of their influence on the retention behaviors of D in irradiated Be.

The initial depth-distribution of D retention in Be with different incident angles, calculated using TRIM code,²⁸ is shown in Fig. 5(a). With the incident angles increasing, the incident D prefers to accumulate near the surface and the total D retention decreases gradually due to the back-scattering effect. The diffusion and reaction effects taken into account in our CD model aggravate the reduction of D retention in Be by further absorbing D near the surface of Be, as shown in Fig. 5(b). The relation of D retention in Be with different incident angles can be fitted by an effective diffusion model as,

where P_t(θ) represents the retention ratio of D in irradiated Be after time t, α is the ratio of D atom diffusion into the bulk (about 0.51), is the average depth of implanted particles which is related to the incident angles, and D_eff represents the effective diffusion coefficient, fitted by the value of ∼9.07 × 10⁻¹⁴ m² s⁻¹. This is reasonable for an effective diffusion coefficient smaller than the value of 1.92 × 10⁻¹³ m² s⁻¹, after including the absorption by different sinks. As shown in Fig. 5(a), the first term of eqn (2) represents the fraction of D diffusing into the bulk, which follows a linear relation fitted to the initial distribution. The second term is the fraction diffusing to the surface, as described by the one-dimensional diffusion theory. From the results, it can be deduced that considering the factor of incident angles, the contribution of the back-scattering effect is greater than that of the long-term diffusion-reaction effect for the D retention in irradiated Be. Additionally, the diffusion-reaction effect can be described by an effective diffusion model, in which the effective diffusion coefficient can be extracted from retention-angle relations. These results should be very helpful for understanding the diffusion and reaction mechanisms during D desorption from bulk Be.

Fig. 5 Concentration profiles of deuterium retention in Be implanted under a flux of 3 × 10¹⁹ D m⁻² s⁻¹ up to a fluence of 1.9 × 10²¹ D m⁻² with an energy of 9 keV at 300 K with the different incident angles, (a) initial depth-distribution of D and (b) considering the long-term dynamical evolution.

Concerning the temperature, the thermal desorption mechanism of D from Be can also be revealed from our results. As shown in Fig. 6, the relative concentrations of D atoms to D₂ molecules are inverse, but the total amount of D retention is almost invariable with increasing temperature until about 400 K, which is in accordance to the experimental fact that the distinct D desorption begins at this temperature (see the inset of Fig. 6).² It can be deduced that the temperature below 400 K is not enough to induce thermal desorption, but improves the reaction rate of D atoms to D₂ molecules. With the temperature higher than 400 K, the amounts of D and D₂ decrease dramatically, due to the high diffusing ability of D and D₂ at higher temperature. In addition, the difference in diffusivity for D atoms (1.92 × 10⁻¹³ m² s⁻¹) and D₂ molecules (1.02 × 10⁻¹⁴ m² s⁻¹) leads to the presence of two thermal desorption peaks at 440 K and 460 K. It is also found that nearly no desorption occurs for immobile D complex clusters below 500 K, due to their high binding energies.

Fig. 6 Detailed temperature profiles of all D, mobile_D defects (D atoms and D₂ molecules), and immobile_D defects in Be under a flux of 3 × 10¹⁹ D m⁻² s⁻¹ up to a fluence of 1.9 × 10²¹ D m⁻² with an incident energy of 9 keV. Inset is the experimental desorption rate with a heating rate of 1 K s⁻¹ under an incident D fluence of 2 × 10²¹ m⁻².²

4 Conclusions

This work has presented the diverse retention behaviors of deuterium in Be under several irradiation conditions by analyzing the competition of diffusion and accumulation. In the high energy incidence case, the forms of D atoms, D₂ molecules and D_mV mainly dominate the D depth-distribution in irradiated Be because of the relatively weak sink strength. The implanted fluence could manipulate the relative proportion of deuterium retention forms between atoms and molecules. The surface desorption, declining incidence and temperature variation also have effects on the behaviors of D retention in irradiated Be and the trapped forms. Additionally, synergistic neutron irradiation induces a uniform distribution of extra point defects which compensates for the recombination zone, but affects the deuterium retention slightly. These results will be helpful in understanding the long-term dynamic evolution mechanics of defects in PFMs.

Acknowledgements

This work was supported by the National Science Foundation of China under Grant No. 11275229, 11475215 & NSAF U1230202, the Special Funds for Major State Basic Research Project of China (973) under Grant No. 2012CB933702, and Director Grants of CASHIPS. Part of the calculations were performed at the Center for Computational Science of CASHIPS, the ScGrid of Supercomputing Center, and the Computer Network Information Center of the Chinese Academy of Sciences.

References

1. D. M. Duffy, Mater. Today, 2009, 12, 38–44.
2. M. Reinelt, A. Allouche, M. Oberkofler and C. Linsmeier, et al., New J. Phys., 2009, 11, 043023.
3. D. M. Duffy, Philos. Trans. R. Soc., A, 2010, 368, 3315–3328.
4. P. B. Zhang, J. J. Zhao and Z. B. Wen, et al., J. Nucl. Mater., 2012, 423, 164–169.
5. V. M. Sharapov, L. E. Gavrilov and V. S. Kulikauskas, et al., J. Nucl. Mater., 1996, 233–237, 870–873.
6. V. M. Sharapov, L. E. Gavrilov and V. K. Alimov, Plasma Devices Oper., 1998, 6, 89–94.
7. H. Krimme and M. Fahnle, J. Nucl. Mater., 1998, 255, 72–74.
8. M. G. Ganchenkova, V. A. Borodin and R. M. Nieminen, et al., Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 134101.
9. M. G. Ganchenkova and V. A. Borodin, Phys. Rev. B, 2007, 75, 054108.
10. P. B. Zhang, J. J. Zhao and B. Wen, J. Phys.: Condens. Matter, 2012, 24, 095004.
11. R. A. Anderl, R. A. Causey and J. W. Davis, et al., J. Nucl. Mater., 1999, 273, 1–26.
12. Y. G. Li, W. H. Zhou, L. F. Huang and Z. Zeng, et al., J. Nucl. Mater., 2012, 431, 26–32.
13. R. H. Ning, Y. G. Li, W. H. Zhou and Z. Zeng, et al., J. Nucl. Mater., 2012, 430, 20–26.
14. T. Faney and B. D. Wirth, Modell. Simul. Mater. Sci. Eng., 2014, 22, 065010.
15. Y. Watanabe, H. Iwakiri and N. Yoshida, et al., Nucl. Instrum. Methods Phys. Res., Sect. B, 2007, 255, 32–36.
16. Y. G. Li, W. H. Zhou and L. F. Huang, et al., Plasma Sci. Technol., 2012, 14(7), 624.
17. R. H. Ning, Y. G. Li, W. H. Zhou and Z. Zeng, et al., Int. J. Mod. Phys. C, 2012, 23, 1250042.
18. Y. G. Li, W. H. Zhou and R. H. Ning, et al., Commun. Comput. Phys., 2012, 11, 1547–1568.
19. C. H. Woo, Philos. Mag. A, 1991, 63, 915–923.
20. F. Christien and A. Barbu, J. Nucl. Mater., 2005, 346, 272–281.
21. F. Christien and A. Barbu, J. Nucl. Mater., 2009, 393, 153–161.
22. T. Faney and S. I. Krasheninnikov, et al., Nucl. Fusion, 2015, 55, 013014.
23. V. N. Chernikov, V. K. Alimov and A. V. Markin, et al., J. Nucl. Mater., 1996, 233–237, 860–864.
24. R. C. Rau, Philos. Mag., 1967, 16, 663–667.
25. C. May and P. Klimanek, II Nuovo Cimento D, 1997, 19, 591–598.
26. S. C. Middleburgh and R. W. Grimes, Acta Mater., 2011, 59, 7095–7103.
27. M. G. Ganchenkova, P. V. Vladimirov and V. A. Borodin, J. Nucl. Mater., 2009, 386–388, 79–81.
28. J. P. Biersack and L. G. Haggmark, Nucl. Instrum. Methods, 1980, 174, 257.
29. N. M. Ghoniem and S. Sharafat, J. Nucl. Mater., 1980, 92, 121–135.
30. M. P. Surh, J. B. Sturgeon and W. G. Wolfer, J. Nucl. Mater., 2004, 325, 44–52.
31. L. Petzold, SIAM J. Sci. Stat. Comput., 1983, 4, 136.
32. B. D. Wirth, G. R. Odette and J. Marian, et al., J. Nucl. Mater., 2004, 329–333, 103–111.
33. C. J. Ortiz, M. J. Caturla and C. C. Fu, et al., Phys. Rev. B., 2007, 75, 100102.
34. V. Chakin, R. Rolli and P. Vladimirov, et al., Phys. Scr., 2009, 138, 014035.
35. C. C. Fu, J. D. Torre and F. Willaime, et al., Nat. Mater., 2005, 4, 68.
36. A. H. Duparc and C. Moingeon, et al., J. Nucl. Mater., 2002, 302, 143–155.
37. V. K. Alimov, V. N. Chernikov and A. P. Zakharov, J. Nucl. Mater., 1997, 241–243, 1047–1052.
38. N. M. Ghoniem, Radiat. Eff. Defects Solids, 1999, 148, 269–318.
39. A. Hu and A. Hassanein, J. Nucl. Mater., 2014, 446, 56–62.

---