Simulation of radiation damages in molybdenum by combining molecular dynamics and OKMC

NUCLEAR CHEMISTRY, RADIOCHEMISTRY, RADIOPHARMACEUTICALS AND NUCLEAR MEDICINE

Simulation of radiation damages in molybdenum by combining molecular dynamics and OKMC

Gui-Yan Wu，

Neng-Wen Hu，

Hui-Qiu Deng，

Shi-Fang Xiao，

Wang-Yu Hu

Nuclear Science and Techniques

Vol.28, No.1

Article number 3

Published in print 01 Jan 2017

Available online 30 Nov 2016

DOI：10.1007/s41365-016-0164-9

34000

In this paper, radiation defects in bcc molybdenum with the PKA (Primary Knock-on Atom) energies of 2–40 keV are simulated by the molecular dynamics. The binding energy of single point defect-to-defect clusters increases with the cluster size. The stability and mobility of point defects and defect clusters are analyzed. The interstitial-type clusters are found to be easily migrating along the <111> direction with low barriers (0.01–0.10 eV). Then, the object kinetic Monte Carlo (OKMC) is used to gain insight into the long-term defect evolution in the cascade. The simulation results indicate that Stage I almost occurs at annealing temperature of 100K, which corresponds to the correlated recombination resulting from the motion of small interstitial clusters (n≤2). The formation of sub-stage partly as result of the small vacancy clusters motion. At about 460 K, the Stage II starts because of uncorrelated recombination due to an emitting mechanism of larger clusters. Size distribution of the clusters at the cascade quenching stage is positively correlated with the PKA energies, affecting notably the subsequent annealing process.

Displacement cascadesMolybdenumDefectsComputer simulation

1. Introduction

With excellent high temperature strength, high melting point, good thermal conductivity, good resistance to corrosion and radiation-induced swelling in nuclear environment, molybdenum (Mo) is a promising material for innovative nuclear energy systems (fast neutron reactors, thorium-fueled molten salt reactors and fusion reactors)^[1-6]. As a structural material of the first-wall^[7,8], Mo will be subject to the impact of high-energy particles, hence the formation of a large number of point defects, and even defects clusters or loops. Further interaction and evolution of the clusters or loops causes microstructure changes and mechanical performance degradation, e.g. irradiation hardening, embrittlement, swelling and creep^[9]. Many features of the radiation damages have been studied by TEM ^[10-13], but further explanations are still inaccessible, such as the mobility of self-interstitials or their clusters, the formation and mobility of the large prismatic vacancy loops and interstitial loops ^[14]. It is, therefore, of significance to understand the radiation damages, including the long time-scale evolution behaviors for defects and defect clusters, and predicting or assessing the material performance in applications.

For defect kinetic evolution in Mo, limited atomistic simulation and experimental data are available for the mobility model of vacancy/interstitial clusters of different sizes. To provide accessible fundamental data for large-scale molecular dynamics simulation of radiation damage, Han et al. ^[15] reported an extensive ab initio study of self-interstitials structures, formation energies, and migration barriers in V and Mo. They found that the <111> dumbbell migration energy in Mo was 0 to within the uncertainty of the current calculation (–0.05 ~ 0.05 eV). This is consistent with experiments (0.05 eV)^[16].

Simulations by Manh et al.^[17] studied the systematic trends for the self-interstitials atom (SIA) defect behavior in bcc transition metal. They found that the <111> configuration was the most stable SIA defect in all the nonmagnetic bcc metals. Their study concludes that group-specific trends showed by the SIA defect structures correlates with the observed thermally activated mobility of SIA defects.

Making clear the formation and evolution of defects and clusters/loops is of practical significance in reducing irradiation embrittlement. Recent studies on displacement cascades produced by 50 keV Xe ions in Mo single crystal revealed that large prismatic vacancy loop formed near the impact surface as a result of fast recrystallization, and the mobility of vacancy dislocation loops was high ^[14].On the contrary, there is no indication that large SIA clusters are formed in Fe in early stage of the cascade evolution ^[18]. A meaningful database of the point defect survival efficiency and partitioning probability into different sized clusters was obtained by Selby et al. ^[19], but mobility of vacancy/SIA clusters in Mo remains unevaluated.

Molecular dynamics (MD) is a useful tool for studying the primary damage produced by high energy displacement cascades and analyzing the migration mechanism of defects. The main characteristics of defect evolution can be predicted using object kinetic Monte Carlo (OKMC). Combining OKMC with MD is a feasible solution to perform atomistic level simulation of the microstructure evolution. This multi-scale technique has been successfully employed to describe the long time evolution of the defects in bcc Fe ^[20,21], fcc Cu ^[22], Ni ^[23], and 3C-SiC ^[24].In their study on evolution of the defect distribution and how the cascade defect distribution evolves in diffusion process, Heinisch and Singh^[25] reported that the spatial segregation, and clustering of the vacancies and SIAs, contributed to a differential production of mobile vacancies and SIAs, with strong temperature dependence. However, the inmost law and mechanisms for defect clusters and loops interactions are still not clear.

In this paper, we report our annealing simulations of radiation-induced damage in bcc molybdenum based on a reliable MD database with PKA (Primary knock-on atom) energies of 2–40 keV. Specifically, we aim to address the following issues: (1) the defect production and the size distribution of defects clusters; (2) the stability and migration behavior of vacancy/SIA clusters; and (3) the annealing of damage as a function of temperature for cascade. This study provides qualitative and quantitative reference information to predict the microstructure evolution of materials under irradiation.

2. Methodology

Molecular dynamics was performed to study displacement cascades in Mo with modified MOLDY codes. A modified analytical embedded-atom method (MAEAM) potential was employed for describing the Mo-Mo interaction. Details of the specific formalisms can be found in Ref. [26].The process of simulating the cascade damage production involved setting up a simulation box of pure molybdenum, sized at 60a₀×60a₀×60a₀ (totally 432,000 atoms), where a₀ is the lattice parameter of body-centered cubic (bcc) Mo bulk. The system was equilibrated at 100K and 300K for 10 pico-seconds. PKA energies of 2, 5, 10, 20, 40 and 60 keV were chosen along the<135>orientation, so as to avoid channeling effect. Periodic boundary conditions were applied for all the MD simulations.

The MD results are essential for OKMC simulations. The migration energies and dissociation energies (the sum of binding and migration energies) of the defect clusters were implemented as key input parameters in kinetic model. The binding energies were obtained by MD method using MAEAM potential ^[27], and the migration energies, by NEB (nudged elastic band) ^[28,29], an efficient method for finding MEPs (minimum energy path) ^[30] between the given initial state and final state (equilibrium configurations). The MEP was found by constructing a set of images (25 images, including two stable images) between the initial and final state.

The OKMC simulation method is capable of tracking the evolution of defects at longer time-scale. In the model, the point defects (vacancies, interstitials and their clusters) were treated as "objects". The objects may migrate, emit and interact with each other, which are defined as "internal event" (migration, recombination, dissociation and annihilation) ^[31-35]. Definitely, are combination event happens when one object encounters its opposite type. If the object is large enough to allow point defect emission, a dissociation event occurs.

The probabilities Γ_n,iof migration jumps and emissions (for object n and event i) are given by ^[36,37]

Γ_{n, i} = ν_{n, i} e x p [- E_{a, n, i} / (k T)],

(1)

where ν_n,i is attempt frequency, E_a,n,i is the corresponding dissociation energy, k is the Boltzmann’s constant and T is the absolute temperature. According to the probabilities Γ_n,i the average time step length is given by ^[38]

Δ t = - l n ξ / k_{t o t},

(2)

where k_tot is the sum of probabilities of all the events, and ξ is a random number in the range between zero and one.

The defect configuration obtained from MD simulations was introduced in the center of the KMC simulation box of Mo bulk, sized at (120a₀)³. The damage annealing simulation was performed for all cascades, and the system was annealed at temperatures from 30 K to 1000 K for a given period of time. The temperature interval was ΔT_j+1=0.05T_j, (j=1,2,3···). The annealing at each temperature was repeated for 50 times. We assumed that larger vacancy clusters (> 4) were immobile.

3, Result and discussion

3.1 The distribution of point-defects

Consistent with previous works ^[19], the number of point defects shortly increased and then decreased due to the interstitials recombining with vacancies. Then, it almost remains stable. It was worth noting that the sub-cascades were observed at the PKA energy of 40 keV. As shown in Table 1, the temperature and PKA energy significantly affect the production and recombination of defects. Higher temperature and PKA energy induced more defects and larger clusters and delayed the time to enter thermal peak stage. At T=100 K, for example, the maximum number of defects is N_d-max=658 and the time to entering thermal peak stage is t_peak =0.342 ps at E_PKA=2 keV; while they are N_d-max=90314 and t_peak =1.662 ps at E_PKA=40 keV. This was probably because the mobility of atoms increases at higher temperature and PKA energy, which hinder recombination of the point defects.

The time to enter thermal peak stage (t_peak), maximun number of defects (N_d-max) and residual Frenkle pairs (N_FP)at different PKA energiesand temperatures.

E_PKA (keV)	t_peak(ps)		N_d-max		N_FP
	100 K	300 K	100 K	300 K	100 K	300K
2	0.342	0.410	658	696	11	15
5	0.408	0.463	1650	1784	48	49
10	0.461	0.840	3434	5167	74	77
20	0.755	1.143	8233	17187	139	199
40	1.662	1.590	90314	60652	567	486

The size (defined by the number of point defects contained) distribution of clusters and defects distribution of the system at T=100K and 300K are shown in Fig.1. It can be seen that most of the clusters were composed of just one or two defects, and large vacancy clusters could be dispersed in the high temperature cascade. Comparing Figs. 1(a–d), the surviving defects clusters generated by high recoil energies are larger than those by the lower energies. In addition, three configurations single interstitials were observed in simulation including<111>oriented crowdions, <110> and <111> direction dumbbells. The most stable interstitial in Mo are the <111> dumbbells, which can easily move along the <111> direction. This is consistent with previous GGA and LDA results ^[14,15].

Fig. 1.

Size distribution of defect clusters at MD temperature of 100 K and 300 K and E_PKA of 10 and 40keV. Blue, vacancy; red, interstitial.

3.2 The migration and binding energies of interstitial/vacancy clusters

The most important parameter for OKMC simulations are the migration energies, which can be obtained by the NEB method, as listed in Table 2. The initial and final states must be specified as searching the saddle point; therefore, the migration direction, equilibrium configuration and migration barrier should be included. In terms of migration direction, the NEB may be only appropriate for single vacancy and small size interstitial clusters. We built the configurations as many as possible, and kept the one having higher binding energy and lower formation energy as the most stable configuration. The migration barriers were calculated by NEB method in a box sized at (10a₀)³, where it was relaxed for 10 ps under 0K. Several equilibrium configurations for vacancy clusters in bcc were defined in Ref. [40, 41]. For example, the migration of mono-vacancy starts from the original point (0, 0, 0) and ended at the nearest neighbor body-centered site (0.5, 0.5, 0.5), by moving along the <111>-direction (Fig. 2a); other clusters were calculated in a similar way. According to Fu et al. ^[39], the migration mechanism of vacancy clusters was simply explained by the fact that the atom at final site migrated to the initial one, recombined with one vacancy and then produced a new vacancy. Furthermore, the motion of small vacancy clusters required passing through metastable intermediate states (for instant, V₂ and V₄, see Figs. 2b and 2d).The same as in other OKMC simulation ^[42], larger vacancy clusters (size>5) were believed to be immobile.

The binding and migration energies of defect clusters.

Defect type	bcc-Mo		α-Fe ^[39] [E_m(eV)]
	E_b(eV)	E_m(eV)
V₁	—	0.985	0.67
V₂	0.268	1.023	0.62
V₃	0.930	0.982	0.35
V₄	1.985	1.383	0.48
V_n	0.055n^2.461	immobile	—
I₁	—	0.009	0.34
I₂	0.349	0.012	0.42
I₃	1.006	0.141	0.43
I₄	1.647	0.113	—
I_n	—	0.0097n^0.0441	—

Fig. 2.

The migration energies of vacancy-type (V₁, V₂, V₃ and V₄) defects moving along the most favorable migration pathways. The x axis is the image points. The vacant sites are represented by green balls.

It should be emphasize that the migration energies of vacancy clusters are size-dependent. For small size V_n clusters, the migration energies are 0.958, 1.032, 0.982 and 1.383eV, respectively, corresponding to n= 1, 2, 3 and 4, as shown in Fig. 2 and Table 2. Unexpectedly, the value of V₃is lower than V₂and V₄ because V₃ cluster has more diversified and complicated configurations. It is interesting that the unexpectedly appearances are similar with ab initio calculations in iron ^[39] (Table 2).One possible reason is that molybdenum and iron are of the same crystal structure. In Table 2, the migration energies of interstitial I_n clusters are much less than that of V_n ones, which illustrates that the migration of I_n clusters is much easier. The migration energies of I_n clusters increased with their sizes, and an empirically-fitted relationship, 0.0097n^0.0441, was fitted. As mentioned before, searching the most stable configurations and exploring several possible migration mechanisms are essential for calculating migration energies. Notably, interstitial-type defects are found to be made of <111> parallel single dumbbells and believed to have special migration mechanisms (Johnson’s mechanism ^[41]), where two atoms of the dumbbells dissociate at first, and then one of them would migrate and extrude the other atom along the migration path to produce a new dumbbell; and the left one would recover at normal lattice site. We found that the I_n clusters could easily migrate with lower migration energies (0.01–0.10 eV). The migration barriers and paths are potted in Fig. 3. All the saddle points lied in middle of the migration path. The binding energies of defect clusters are the other energetic input data for the simulations. Their values (V_2-4 and I_2-4) are determined by MD simulations with MAEAM potential. Table 2 shows that the binding energies increase with the size. The stable structures are shown in Fig. 2.

Fig. 3.

The migration energies of interstitial-type (I, I₂, I₃ and I₄) defects along the most favorable migration pathways. The x axis is the image points.

3.3 The OKMC simulation

In the OKMC simulation, the displacement cascades simulated by MD (Table 2) were used as the initial defect distribution. Cascade energies of 2, 5, 10, 20 and 40 keV were annealed. Each cascade result was introduced in the center of the OKMC box sized at (120 a₀)³. During annealing, the system was simulated at temperatures of 30–1000K and each cascade was annealed 50 times. Fig. 4 shows the average number of defects as the function of annealing temperature. It is interesting that each annealing curve has typical stages for the decrease in the number of defects Stage I is around 100 K. This is the recombination stage, mainly due to the motion of interstitials clusters which were annihilated with nearby vacancy defects or larger clusters. Stage II, at 390–460K, is related to the recombination of the uncorrelated defects. For 5 keV and 10 keV cascades formed at 100 K, Stage II occurred at 390–395 K; while it occurred at 430 K for the100 K 40 keV cascade. There exist two mechanisms acting together on the uncorrelated recombination; the irrelative recombination of small size defects (<3), and the defects emitting from the larger size clusters. The 40 keV results are mainly due to the latter.

Fig. 4.

Residual number of vacancies vs. annealing temperature, with primary defects formed at T=100 K or 300 K and E_PKA=10 keV or 40 keV.

The most significant feature is a slight fluctuation observed in 300–600 K (Fig. 5a), where the number of defects increases first but then decreases. It declines continuously and finally reaches a stable one. This is mainly due to the emitting events during the annealing process, in which a large number of new tiny clusters are extracted, and most of them migrate to annihilate with opposite defects. Fig. 5 also shows that for 2 keV cascade, the second recombination of the defects started at 537K, being100 K higher than those of the 20keV cascade. This is because that the small clusters generated at E_PKA= 2keV completely recombine in Stage I, but larger clusters recombine at high temperatures. Note that most of small clusters occur at higher E_PKA, especially in Stage II. So, the temperature reduced by the energy augments.

Fig. 5.

Residual number of vacancies vs. the annealing temperature, for primary defects at 300 K and E_PKA=2 and20 keV.

As shown in Fig. (1), the defect distribution and cluster size are affected remarkably by PKA energies, and so is the recombination efficiency in annealing. Fig. 6 shows the fraction of recombination as a function of annealing temperature. Apparently, a higher cascade temperature leads to a higher ratio of recombination, but the impact of cascade energy is otherwise. This indicates that the cascade energy is notably positively correlated with the cluster size and thus affects the annealing process. For 80% defects formed at E_PKA =10 keV are recombined before the annealing temperature goes to around 100K, due to mainly the small size interstitial clusters migration along <111> direction. In particular, Stage I has several sub-stages, because that as the annealing temperature increases, more small size vacancy clusters move and recombine with mobile interstitial clusters. Above 460K, the recovery ratio fluctuates slightly because almost all small interstitial clusters are recombined and fewer mobile vacancy clusters breakaway from larger clusters. In Stage II, an additional 30% recombination ratio are added for cascade 2 keV and 10 keV at 100K (cascade temperature). This is mainly due to the motion of residual small vacancy clusters and the emitting events during the annealing process. From Fig. 6, the fractions of surviving defects are20%–75% as Stage I ends, and 0–70% as the annealing ends. So, once again, the size and number distribution of clusters from the MD simulation can make a great difference in annealing simulation.

Fig. 6.

The fraction of recovery defects as function of annealing temperature for E_PKA = 2, 10 and 40 keV.

Figure 7 shows the average number of clusters for a 10 keV cascade, as the function of annealing temperature. It can be seen that a rapid drop of the number of small size interstitial clusters (n≤2) at the beginning of Stage I, and the number of large size interstitial clusters (2<n≤10) plunged. At this stage, large size vacancy clusters can be annihilated with the fast-moving interstitials clusters (n≤2), and the shrunken vacancy clusters further annihilate with large size interstitial clusters (2<n≤10). Therefore, the results illustrate that the correlated recombination stage resulting from the motion of small interstitial clusters to a great extent. Also, the number of small vacancy clusters (n≤2) decreases, hence the occurrence of sub-stages. At460 K, at the uncorrelated recombination stage, the number of small size vacancy clusters (n≤2) and the large size interstitial clusters (2<n≤10) begin to predominantly decline. This is due to the enhanced motion of small size vacancy clusters and increased possibility of emitting event. In addition, from the MD simulations, the interstitials prefer migrating along the <111> direction, but the high annealing temperature may change their preference and they may recombine with the vacancy clusters in the core.

Fig. 7.

The average number of species clusters versus annealing temperature for PKA energy of 10 keV.

4. Conclusion

MD and OKMC simulations were applied to investigate the displacement cascades and the evolution of irradiation-induced defects in Mo. The number of defect surviving fraction was 0.02–0.05, indicating that most point defects recombined before annealing. Vacancy clusters were observed and the <111> direction dumbbells were proved to be the most stable interstitials in Mo and could move along the <111> direction with low migration barriers (0.01–0.10 eV). Stage I at about 100 K was attributed to the correlated recombination, results from the motion SIAs (n≤2) to a great extent. The possibility of sub-stage formation as due to the small vacancy clusters motion. In addition, a slight fluctuation was observed at300–600 K. This was mainly due to the emitting events during the annealing. The size distribution of the clusters at the cascade quenching stage was positively correlated with the PKA energies, affecting notably the subsequent annealing process.

References:

[1]

P Yvon and F Carré.

Structural materials challenges for advanced reactor systems

. J. Nucl. Mater, 2009, 385: 217-222.DOI: 10.1016/j.jnucmat.2008.11.026