Radiation-induced damage and evolution of defects in Mo

PHYSICAL REVIEW B 84, 104109 (2011) Radiation-induced damage and evolution of defects in Mo Sergey V. Starikov,* Zeke Insepov,† and Jeffrey Rest Argonne National Laboratory, Argonne, Illinois 60439, USA Alexey Yu. Kuksin, Genri E. Norman, Vladimir V. Stegailov, and Alexey V. Yanilkin Joint Institute for High Temperatures, Russian Academy of Sciences, Moscow 125412, Russia and Moscow Institute of Physics and Technology, Dolgoprudny 141700, Russia (Received 30 December 2010; revised manuscript received 13 May 2011; published 8 September 2011) The formation of defects in bcc Mo lattice as a result of 50-keV Xe bombardment is studied via atomistic simulation with an interatomic potential developed using the force-matching ab initio based approach. The defect evolution in the cascade is described. Diffusion and interaction of interstitials and vacancies are analyzed. Only small interstitial atom clusters form directly in the cascade. Larger clusters grow only via aggregation at temperatures up to 2000 K. Stable forms of clusters demonstrate one-dimensional diffusion with a very high diffusion coefficient and escape quickly to the open surface. Point vacancies have much lower diffusivity and do not aggregate. The possibility of a large prismatic vacancy loop formation near the impact surface as a result of fast recrystallization is revealed. The mobility of the vacancy dislocation loop segments is high, however, the motion of the entire loops is strongly hindered by neighbor point defects. This paper explains the existence of the large prismatic vacancy loops and the absence of the interstitial loops in the recent experiments with ion irradiation of Mo foils. DOI: 10.1103/PhysRevB.84.104109 PACS number(s): 61.80.−x, 64.30.−t, 64.70.D−, 71.10.Li I. INTRODUCTION The strength and mechanical properties of nuclear fuel pellets as well as structural materials of nuclear reactors essentially depend on the influence of neutron irradiation and the corresponding evolution of the material microstructure. The evolution of radiation defects is experimentally studied using transmission electron microscopy in the samples subjected to neutron irradiation or heavy-ion bombardment.1–8 The method of heavy-ion implantation mimics the primary knock-on atoms (PKA) produced by collisions with neutrons. This method provides more control of the irradiation conditions. However, the penetration depth of ions is smaller in comparison with neutrons, and the ion cascades start near the ion-entry surface (the influence of surfaces is especially pronounced in experiments with thin foils, e.g., see Ref. 8). The detailed knowledge of the ion cascade evolution is crucial for the application of the heavy-ion implantation data for neutron irradiation conditions. Although many features of the radiation damage have been explained, several phenomena are still discussed, including the mechanisms of production of large clusters of defects and dislocation loops,5,8–12 the effect of PKA mass on the defect yield,8,12 the influence of surfaces in the experiments with thin foils,8 the mobility of self-interstitial atoms (SIAs), clusters (SIACs), and dislocation loops.2,4,13–16 One of the most plausible explanations of the formation of superlattices of bubbles1–3,13,14,16 in several metals under the radiation damage is based on the anisotropic diffusion of SIAs or SIACs. However, there are contradicting statements in the previous works about the mechanisms and dimensionality of di-interstitial SIA (di-SIA) and SIAC diffusion.2,13,14,16,17 At the same time, the comparison with the recent experimental results showed a large discrepancy in the diffusivity of dislocation loops.4 1098-0121/2011/84(10)/104109(8) The understanding of the mechanisms of the production of large clusters of defects is of practical importance for reducing the irradiation embrittlement. The presence of impurity atoms and defect sinks can substantially influence the processes of the radiation softening or hardening if the sessile clusters are formed via diffusion-controlled nucleation and growth that has been confirmed for some metals (Ni, Cu, Ag, Pt, and W). However, the formation mechanism of SIACs in Al, Fe, and Mo remains still unclear. In the case of the in-cascade formation of clusters, the influence of impurities is not clear.18 Recent experiments on low-dose (<1016 ions/m2 ) irradiation of Mo single-crystal foils revealed the formation of vacancy loops within individual cascades by a process of cascade collapse.8 In iron, under self-ion irradiation, the cascade collapse is not observed at low doses,6 while visible damage is seen at high ion doses when cascade overlap becomes significant as well as in the cascades initiated by more massive ions.7 For both Mo and Fe, it was shown that the number of visible loops increases as the mass of the incident ion becomes larger,6–8 e.g., in the case of molecular ion irradiation. Dislocation loops mostly composed of vacancies are observed. Interstitial loops have high mobility, and they are believed to go to the surfaces and disappear. A useful tool for studying damage of materials at high-energy deposition,15,19–23 including displacement cascades,10–12,18,24–30 is a classical molecular dynamics (MD) method. A large amount of the displacement cascades simulations for bcc metals was done on the example of Fe.10–12,18,24,25,27–29 Cascades in Mo were studied in Ref. 26. Most of the MD results on radiation damage were obtained within the PKA model, and the influence of the surface attracted less attention.30 Both classical MD and electron density functional theory (DFT) atomistic models were applied to analyze the migration 104109-1 ©2011 American Physical Society PHYSICAL REVIEW B 84, 104109 (2011) SERGEY V. STARIKOV et al. mechanisms of defects and to evaluate their mobility (e.g., see Refs. 15, 31, and 32,). The resulting data can be deployed as an input for the theoretical models of long time evolution of damaged material.3,14,33,34 It was shown that both the SIAs and the vacancy loops have high one-dimensional (1D) mobility along close-packed crystal directions in bcc crystals.10,31,32 The small clusters can switch their propagation direction between equivalent directions, which result in the three-dimensional (3D) diffusion at elevated temperatures. While the influence of the interatomic potential in classical MD cascade models is still not quite clear,11,25 it is concluded that the number and sizes of vacancy clusters produced during cascade collapse in Fe under self-ion irradiation are small in a qualitative agreement with experimental observations. The increase of the PKA mass results in formation of large SIA and vacancy clusters or loops,12,29 which are potentially resolvable by the transmission electron microscopy. However, a certain disagreement exists. Most loops observed in the thin-foil experiments are composed of vacancies.7 This peculiarity is more pronounced for Mo.8 It can be only partially attributed to the fast glide of SIA loops to the open surfaces8 since the foils with the 110 surface should preserve a certain fraction of SIA loops with a glide cylinder aligned in the foil plane. The defect evolution described above is studied in this work for the molybdenum single crystal. Mo is a refractory metal with high strength, radiation, and creep resistance at high temperatures. The study of radiation damage in pure Mo can shed light on the mechanical properties of nuclear fuel based on U-Mo alloys. Basic properties (such as formation energies and structure) of SIAs in Mo were studied using DFT.31,32 It was found that the difference in formation energies of 111 crowdion and dumbbell configurations is negligibly small, which implies the high SIA mobility along the 111 direction at low temperatures in agreement with the resistivity recovery experiments.35 The old generation interatomic potentials for Mo26,36 do not reproduce correctly the hierarchy of formation energies of SIAs that results in the absence of 1D diffusion and very small SIA diffusivity at low temperatures. The potential proposed in Ref. 37 takes into account this hierarchy but predicts anomalous thermal expansion and gives a large overestimation of the melting temperature. In this paper, we construct an interatomic potential for Mo that correctly describes thermodynamic properties and energies of the main defect structures and apply this potential to the radiation cascade studies. The interatomic potential and the procedure of its construction are described in the Sec. II. Section III reviews the details of the cascade simulations. Section IV is devoted to the processes of defect interaction. In Sec. V, the formation of vacancy loops near the open surface is discussed. II. INTERATOMIC POTENTIAL Here, we develop an embedded-atom method38 (EAM) interatomic potential for the Mo-Xe system,39 which describes energetic properties of point defects and defect clusters in consistency with DFT calculations. The energy Ui of an atom i is given by   Ui = Fα (ρ i ) + ϕ αβ (rij ),ρi = ρβ (rij ), (1) j =i j =i where F is the embedding energy that is a function of the electron density ρ, ϕ is a pair potential interaction, and α and β are the element types of atoms i and j . The multibody nature of the EAM potential is a result of the embedding energy term. Both summations in the formula are over all neighbors j of atom i within the cutoff distance. All the functions are represented as cubic splines with 10 nodes each. The force-matching method40,41 is used for the development of the potential as implemented in the POTFIT code.41 This method provides a way to construct physically justified interatomic potentials excluding the experimental data from the fitting database. The idea is to adjust the interatomic potential functions to optimally reproduce per-atom forces (together with total energies and stresses) computed at the ab initio level for a fine-tuned set of reference structures. The reference data are calculated using the DFT code VASP (Ref. 42) with approximatly 200 atoms in the supercell (the number of atoms depends on the strcture, the number of defects, etc.) The projector augmented-wave (PAW) pseudopotentials together with the plane-wave basis cutoff energy of 400 eV are used. The Brillouin zone is sampled with the 2 × 2 × 2 Monkhorst-Pack k-point mesh. The generalizedgradient approximation (GGA) of Perdew and Wang43 is used for the exchange-correlation functional. We use 81 configurations with 10 746 atoms all together. These configurations represent 39 states with pure Mo (liquid and bcc solid at different densities; lattice with SIAs and/or FIG. 1. (Color online) The functions of the EAM potential for the Mo-Xe system. Pair potentials: (a) solid red line, the Mo-Mo interaction; black dashed-dotted line, the Mo-Xe interaction; blue short dashed line, the Xe-Xe interaction. Electron densities (b) and embedding functions (c): solid red line, Mo; blue dashed line, Xe. 104109-2 PHYSICAL REVIEW B 84, 104109 (2011) RADIATION-INDUCED DAMAGE AND EVOLUTION OF . . . TABLE I. The comparison of EAM interatomic potentials of molybdenum at room temperature using the cohesive energy Ec , the equilibrium lattice parameter a, the elastic constants Cij , and the melting temperature Tm . The experimental data are taken from Refs. 46, 47, and 52. Ec (eV) a (Å) C11 (GPa) C12 (GPa) Tm (K) Experiment EAM (Ref. 36) EAM (Ref. 37) EAM (this work) 6.82 3.1472 465 176 2890 6.85 3.1471 420 170 6.84 3.1465 600 169 3350 6.88 3.1470 550 220 2630 vacancies and/or surface), 20 states with pure Xe (liquid and solid), and 22 states with Mo and Xe (including a single Xe atom in pure Mo). The energies and stresses for different configurations of point defects are included into the target function with higher weights. Initially, all the atomic configurations are taken from classical MD runs at different temperatures performed with some preliminary EAM potential. The fitting procedure consists of the following steps: (1) fitting the parameters of the potential to the ab initio database, (2) testing the potential with respect to certain properties, and (3) recalculating the initial set of configurations with the updated potential. This procedure is performed in an iterative manner in order to improve the description of the desired properties. The developed potential reproduces energies, stresses, and per-atom forces of the reference data with the average accuracy 0.5%, 4%, and 25%, respectively. The potential functions are shown in Fig. 1. For simulations of high-energy collisions, the pair parts of the EAM potential are modified at small distances to transform smoothly into the universal screened Coulomb potential.44 The basic properties of pure molybdenum calculated with this and two previously published EAM potentials are summarized in Table I and compared with the experimental data. The accuracy of the potential is slightly lower than that of its predecessors at room temperature and zero pressure. However, this potential has significant advantage in the description of processes in a wide range of pressures and temperatures. This EAM potential reproduces the room-temperature isotherms of Mo and Xe up to approximately 600 GPa in agreement with the experimental data.45,46 As shown in Fig. 2, the developed potential reproduces the isothermal compression of pure Mo with good precision. Thermal expansion at zero pressure is reproduced well47 up to the melting point with 1%–2% accuracy as shown in Fig. 3. Such a behavior is required when temperature dependencies of certain physical properties are studied. The melting of pure molybdenum is investigated using twophase simulation.19,48,49 The calculated melting line of pure Mo is very close to the DFT results of Ref. 50 and agrees well with the shock-wave test data.51 Figure 4 shows the melting temperatures Tm at various pressures P . There is a considerable scatter in the experimental results obtained using the diamond anvil cell and shock-wave techniques. The possible reasons for this discrepancy for molybdenum and other metals were discussed in Refs. 19, 50, and 52–54. Relative energies of point-defect configurations and clusters are important for the proper description of their diffusion mechanisms. One can judge about the stability of defects by their formation energy FIG. 2. (Color online) The isotherm of Mo at room temperature: 1, experimental data (Ref. 46); 2, 3 and 4, EAM data for Refs. 37 and 36 and the potential developed in this work, respectively. FIG. 3. (Color online) The thermal expansion of Mo at zero pressure: 1, experimental data (Ref. 47); 2, 3 and 4, EAM data for Refs. 37 and 36 and the potential developed in this work, respectively. Ef = Edef (N) − N Eid (Nid ). Nid (2) Here, Edef is the energy of the system with defect comprising N Mo atoms, and Eid is the total energy of an ideal lattice with Nid atoms. In order to minimize the energy of the system, the relaxation of the atomic positions is performed via the conjugate gradient algorithm. The size of the box is large 104109-3 PHYSICAL REVIEW B 84, 104109 (2011) SERGEY V. STARIKOV et al. FIG. 4. (Color online) The melting line of Mo: 1, the melting points obtained in diamond anvil cells (Ref. 52); 2, the melting point obtained in shock-wave experiment (Ref. 51); 3 and 4, the calculated melting lines from Refs. 50 and 53; 5, two-phase MD simulation with the potential developed in this work. enough (Nid = 250) and the elastic energy due to the defect is negligible. This potential gives the following formation energies of Mo SIAs in different configurations: crowdion 111, 6.42 eV (C111 ); dumbbell 111, 6.43 eV (D111 ); dumbbell 110, 6.67 eV (D110 ); tetrahedral, 7.19 eV (T ); dumbbell 100, 7.35 eV (D100 ); and octahedral, 7.77 eV (O). The vacancy formation energy is 2.79 eV. The relative formation energies of SIA (with respect to a minimum value) are shown in Fig. 5 in comparison with the data from DFT and other EAM potentials. The most stable configurations of SIA are a dumbbell 111 and a crowdion 111 with very small difference in formation energies, in accordance with DFT calculations.31,32 Such a small difference between configurations results in a 1D diffusion of SIA at very low temperatures in agreement with the resistivity recovery measurements following electron irradiation.35 The 1D diffusion proceeds via sequential transformations of SIA from crowdion to dumbbell configurations along the 111 axis. The corresponding SIA diffusivity is high: ∼10−4 cm2 /s at temperatures T  300 K. FIG. 5. (Color online) The relative formation energy of SIA (relative to the minimum value): 1 and 2, DFT calculations in Refs. 31 and 32; 3, 4, and 5, EAM data (minimization) for the potentials from Refs. 37 and 36 and for the potential developed in this work, respectively. FIG. 6. (Color online) The types of di-SIA in Mo: (a) stable mobile type; (b), (c) instable immobile types. The used signs: black solid squares, atoms of the lattice; red circles, interstitial atoms; open square, a missing atom of the lattice. As temperature increases, the 111-110 dumbbell transitions are activated providing rotations of crowdion propagation directions. The height of the corresponding barrier can be estimated as a difference in formation energies of the 111 and 110 dumbbells. Such events are quite frequent on the nanosecond time scale already at 500 K. The dumbbell pair can change its orientation to any of the four equivalent 111 axes. It provides a viable mechanism for 3D SIA diffusion. The average length of 1D segments (the distance between subsequent reorientations) approaches interatomic distances at T ∼ 2000 K. This potential is tested by the reproduction of formation energies of di-SIAs in Mo. Atomic structures of the basic types of di-SIA are shown in Fig. 6. The formation energies of di-SIAs are computed by energy minimization in DFT and MD calculations with this potential. N = 252 atoms are used in both types of calculations. The results are shown in Table II. This potential reproduces correctly the di-SIA energy hierarchy. The more stable type (a) is characterized by 1D diffusion similar to a single SIA. Note that the metastable types of di-SIA (b and c) may exist at low temperatures because of the nonzero energy barrier between them. The formation energies of the Xe inclusion at different positions in Mo lattice are shown in Table III. The lowestenergy type of atom arrangement is the substitutional Xe. It is even favorable for a Xe atom to replace a host Mo atom, thus producing a SIA. A transition of a Xe atom to an interstitial position (either octahedral or tetrahedral) requires considerable energy for activation. This potential may be useful for the simulation of Xe atoms evolution in the Mo matrix (e.g., the bubble formation process) as well as for the simulation in pure Mo. In this work, we study the evolution of Mo structure under high-energy Xe bombardment. III. DETAILS OF CASCADE SIMULATIONS The evolution of radiation damage in displacement cascades generated by a single Xe ion in a Mo lattice is studied using MD simulations. The ion energy is taken to be 48.5 keV. It is assumed that there is no difference between the surface 104109-4 PHYSICAL REVIEW B 84, 104109 (2011) RADIATION-INDUCED DAMAGE AND EVOLUTION OF . . . TABLE II. The calculated formation energies of di-SIAs (the absolute value and the relative value with respect to the SIA doubled energy). The types of di-SIAs are shown in Fig. 6. Ef (DFT) (eV) Ef (EAM) (eV) Type a Type b Type c Doubled energy of SIA 12.42 (−1.36) 11.71 (−1.12) 12.90 (−0.88) 11.92 (−0.91) 13.53 (−0.25) 12.46 (−0.37) 13.78 12.83 bombardment by a Xe atom or a Xe+ ion at such high kinetic energy. The x axis is normal to the impact surface. Periodic boundary conditions are used only in the y and z directions. Two different crystal orientations of the impact surface are considered, {100} and {110}, with the latter chosen in order to provide some of the equivalent 111 glide directions parallel to a free surface. The ion impact direction is about 6◦ offnormal to the open surface. The simulation box varies in size from 30 to 80 lattice constants a in the y and z directions and from 80a to 300a along the x axis. The initial temperature of the target is 250–2500 K. 50 cascades are calculated and analyzed in total. MD runs are carried out by the LAMMPS code.55 The time step is adjusted during the simulation in order to fulfill the criterion of the maximum allowed displacement of atom. During the penetration of a Xe ion into the bulk crystal, the generation of several subcascades takes place. The penetration depth varies significantly from one run to another with the average value of 15 nm. The energy deposited in the subcascades leads to local disordering within several tenths of picoseconds followed by recrystallization. This cascade collapse leads to the formation of vacancies close to the initiating Xe ion trajectory and SIAs on the periphery. A single Xe ion cascade produces ∼102 SIAs and vacancies with smaller numbers at higher T . However, most defects recombine after 1–2 ps. Single vacancies are practically immobile in all the temperature ranges considered (even on the ns time scale). 0.5–1.0 ns after the SIA vacancies recombination, the interaction of moving defects leads to the formation of mobile and sessile SIACs (Figs. 7, 8, and 9). IV. DEFECT INTERACTION IN THE CASCADES We distinguish the following features of the defect interaction process: (i) The SIACs produced directly in the cascade are quite small (2–4 SIAs). A possibility of formation of larger SIA clusters was reported in previous cascade simulations in Mo,26 however, their average size is also quite small. The growth of larger clusters (up to 6–7 SIAs) takes place on the nanosecond time scale as a result of SIA diffusion toward each other. Therefore, large SIACs in Mo are produced only via diffusion and coalescence of SIAs. It is consistent with the experimental evidences obtained for neutron-irradiated Mo,5 where the influence of the surfaces is small and formation of vacancy loops on the surfaces (as discussed below) is neglible. (ii) A SIAC is more stable than separate SIAs. Static calculations show that the bonding energy of the cluster increases with cluster size (like in iron10 ). Therefore, at the collision of two SIAs (or SIA and SIAC), aggregation takes place. The clusters are formed up to T ∼ 2000 K. (iii) Stable forms of SIACs show 1D thermal diffusion (the corresponding diffusion coefficient is of the same order of magnitude as for single SIAs). The most stable form of SIAC corresponds to a bundle of crowdions aligned along the 111 direction see Figs. 7(b), 8(c), and Ref. 56, thus retaining 1D migration (see Ref. 57). Such clusters can be regarded as nuclei of 111 prismatic dislocation loops. The 1D diffusion of dislocation loops has been experimentally observed in iron.4 (iv) At T  1000 K, there exists the possibility of formation of stable immobile SIACs [Figs. 8(a) and 8(b)]. However, they are less stable than mobile SIACs of the same size, and relaxation to mobile structures [Figs. 8(b) and 8(c)] is observed on the ps or ns time scale depending on temperature. This fact can be explained by the existence of an energy barrier between mobile and immobile forms of SIACs. At T  1000 K, only a small amount of immobile SIACs are present. The only possibility for the formation of immobile SIACs at high temperatures is from collisions of large mobile SIACs. The presence of immobile defect clusters at low temperatures agrees with the results for iron.15 (v) Some part of mobile SIAs and clusters migrate to the crystal surface and leave the lattice. An excess concentration of vacancies arises on the ns time scale. In general, the formation of mobile or immobile SIACs from SIAs is equally probable, but at T  1000 K, an immobile SIAC quickly relaxes to mobile structure and can leave the system if its propagation direction is favorable. The number of lost defects is several times smaller for the {110} than for {100} open surface. Therefore, an open surface is an active sink for SIAs and SIACs at high T . V. VACANCY LOOPS FORMATION The formation of a large prismatic vacancy loop is observed in some runs (approximately in 3 cases out of 50 considered) TABLE III. The formation energies for the interstitial and substitutional Xe atom in the bcc Mo lattice. Ef (DFT) (eV) Ef (EAM) (eV) Octahedral Tetrahedral Substitutional Xe + SIA Substitutional Xe 17.1 17.1 16.2 16.5 15.8 15.3 8.8 8.4 104109-5 PHYSICAL REVIEW B 84, 104109 (2011) SERGEY V. STARIKOV et al. FIG. 8. (Color online) Examples of SIACs: (a) an immobile cluster of 15 SIA; (b) a cluster of 13 SIA initially immobile and (c) mobile after reorientation of 3 SIAs. (These structures are specially prepared and not taken from the MD cascade simulations.) FIG. 7. (Color online) The evolution of a displacement cascade in Mo produced by a Xe ion. Atoms comprising defects (blue open circles for vacancies, red solid circles for SIAs) in one half of the simulation box are shown in xy projection at 5 ps and 3 ns after the Xe surface impact (T = 1200 K). The inset shows the main SIA types: (a) a crowdion, (b) 111, and (c) 110 dumbbells. The arrow shows the path of the Xe ion. where the displacement cascade produced remains close to the frontal surface (or sometimes even the rare surface in case of long cascades) (see Fig. 10). The ion stops in the vicinity of the crystal surface and its energy is deposited into a relatively small volume. An intensive expansion of the melt into vacuum decreases the density of the melt region. Compression of the disordered region takes place when its temperature decreases. The recrystallization starts from the bulk at the same time. A few crystalline layers of atoms can form a hill on the surface that leads to a large amount of vacancies left in the bulk. Formation of a prismatic vacancy loop is possible in such a vacancy-rich region at the last stage of recrystallization.58 Moreover, a number of single vacancies can remain around this dislocation loop. Similar results were reported in literature for the Ni3 Al alloy59,60 and discussed hypothetically as an interpretation of the influence of surfaces on defects production in experiments.8 Dislocation loops with Burgers vectors 111 are observed for surfaces of both orientations considered. The image force due to the open surface results in a loop glide along the 111 direction towards the {100} open surface. Near the {110} open surface, the vacancy loop stays trapped by the surrounding point defects within (at least) several hundreds of picoseconds, while some of its segments exhibit 1D oscillations. Such an influence of point defects may be the reason for the difference in the estimates of thermal diffusivity obtained for the dislocation loops in MD and in experiments.4 One should note a significant difference of the presented results with a recent MD simulation12 of cascades produced by Bi ions irradiation of iron (both atomic and molecular): the formation of large SIA clusters was reported in Ref. 12, which is not observed in our simulations. The reason is probably connected with the difference in masses of colliding particles: a Bi atom is much more massive than a Fe atom, while the difference between Xe and Mo atoms is smaller. The ratio of masses in the case of Xe irradiation of the Mo target is very close to one analyzed in experiments8 where Mo is irradiated with Sb ions. The absence of large SIA clusters reported here agrees well with the experimental observations: only vacancy loops were reported in Ref. 8. Such observations can be explained if the vacancy loops are produced in near-surface cascades via the mechanism described above. Then, the enhanced production of the loops in case of molecular ion irradiation (at the same energy per ion) can be connected to the smaller penetration depth of molecular ions in comparison to a single ion: molecular ion cascades form closer to the surface and there is an increased probability of loops formation in a near-surface vacancy-rich region after fast recrystallization. FIG. 9. (Color online) The time dependence of interstitial defect concentrations with respect to the initial vacancy concentration for temperatures (a) 250 K, (b) 800, and (c) 1100 K: 1, SIA; 2, di-interstitial (immobile); 3, di-interstitial (mobile); 4 and 5, clusters with more than 2 SIAs (immobile and mobile). 104109-6 PHYSICAL REVIEW B 84, 104109 (2011) RADIATION-INDUCED DAMAGE AND EVOLUTION OF . . . FIG. 10. (Color online) Vacancy loop formation at radiation damage: (a) t = 5.5 ps; (b) t = 89.3 ps. The solid line is the position of free surface. Only the disordered atoms (atoms forming defects) are shown. T = 1100 K. The dark blue (dark gray) particles denote atoms with low coordination number. The blue/gray and yellow/light gray particles correspond to high coordination number. The other aspect that favors the survival of vacancy loops is that they are formed in a region with large concentration of point defects, which leads to a reduced mobility of such loops. Therefore, the comparison of the current MD simulations with the experimental data on Mo (Ref. 8) shows that the effects of the surface are significant in thin-foil experiments at ion energies ∼102 keV due to the possibility of near-surface production of large vacancy loops. cascades produced by Xe ions with energies about 50 keV in Mo single crystal with open surfaces are studied. Interstitial clusters produced directly in the cascade are quite small (2–3 SIAs). No evidence of the large SIA cluster formations in the early stage of cascade evolution is found, contrary to the observations in Fe.12 The growth of clusters is possible on the nanosecond time scale as a result of SIA diffusion toward each other; however, most of the SIAs possess fast 1D diffusion and can be easily lost on the surfaces. If the energy density deposited by the ion in the vicinity of the open surface is high enough for melting, then the spontaneous nucleation of large vacancy dislocation loop is possible as a result of fast recrystallization. Thus, the nature of the dislocation loops observed agrees well with the experiments8 on ion irradiation of Mo foils. Enhanced production of the loops for molecular ions (at the same energy per ion) is connected with the smaller penetration depth of molecular ions in comparison to single ions. The mobility of the segments of the vacancy dislocation loop is high, however, the motion of the entire loops is strongly hindered by a number of point defects in the region where they have been produced. Thus, it leads to the reduced mobility of such loops and their preferential survival. ACKNOWLEDGMENTS The interatomic potential for Mo is presented, which takes into account the properties of the main defect structures of SIAs and vacancies and reproduces thermodynamic properties in a good agreement with experimental data. Displacement Simulations were carried out on the computing cluster Fusion and the IBM Blue Gene/P at Argonne National Laboratory and on the MIPT-60 cluster of the Moscow Institute of Physics and Technology. This work was supported in part by the US Department of Energy Office of Advanced Scientific Computing Research, Office of Science, under Contract No. DE-AC02-06CH11357, by the Program for Basic Research of the RAS No. 2, Grant No. RFBR 09-08-01116, and the President RF Grants No. MK-3174.2011.8 (A.Yu.K.) and No. MK-6719.2010.8 (V.V.S.). * 12 VI. CONCLUSIONS Also at Joint Institute for High Temperatures of Russian Academy of Sciences and Moscow Institute of Physics and Technology; starikov@ihed.ras.ru. † insepov@anl.gov 1 J. H. Evans, Nature (London) 229, 403 (1971). 2 J. H. Evans, Philos. Mag. Lett. 87, 575 (2007). 3 J. Rest and G. L. Hofman, J. Nucl. Mater. 277, 231 (2000). 4 K. Arakawa, K. Ono, M. Isshiki, K. Mimura, M. Uchikoshi, and H. Mori, Science 318, 956 (2007). 5 M. Li, N. Hashimoto, T. Byun, L. Snead, and S. Zinkle, J. Nucl. Mater. 367-370, 817 (2007). 6 C. A. English, M. L. Jenkins, and B. L. Eyre, Philos. Mag. 38, 97 (1978). 7 Z. Yao, M. Hernandez-Mayoral, M. L. Jenkins, and M. A. Kirk, Philos. Mag. 88, 2851 (2008). 8 C. English and M. Jenkins, Philos. Mag. 90, 821 (2010). 9 T. Diaz de la Rubia and M. W. Guinan, Phys. Rev. Lett. 66, 2766 (1991). 10 N. Soneda and T. D. de la Rubia, Philos. Mag. A 78, 995 (1998). 11 L. Malerba, J. Nucl. Mater. 351, 28 (2006). A. F. Calder, D. J. Bacon, A. V. Barashev, and Y. N. Osetsky, Philos. Mag. 90, 863 (2010). 13 H. L. Heinish and B. N. Singh, J. Nucl. Mater. 307, 876 (2002). 14 A. Semenov, C. Woo, and W. Frank, Appl. Phys. A 93, 365 (2008). 15 D. A. Terentyev, T. P. C. Klaver, P. Olsson, M.-C. Marinica, F. Willaime, C. Domain, and L. Malerba, Phys. Rev. Lett. 100, 145503 (2008). 16 A. V. Barashev and S. I. Golubov, Philos. Mag. 90, 1787 (2010). 17 J. H. Evans, Philos. Mag. 86, 173 (2006). 18 G. Lucas and R. Schaublin, J. Phys. Condens. Matter 20, 415206 (2008). 19 S. V. Starikov and V. V. Stegailov, Phys. Rev. B 80, 220104 (2009). 20 A. Y. Kuksin, G. E. Norman, V. V. Pisarev, V. V. Stegailov, and A. V. Yanilkin, Phys. Rev. B 82, 174101 (2010). 21 D. K. Belashchenko, High Temp. 48, 646 (2010). 22 B. J. Demaske, V. V. Zhakhovsky, N. A. Inogamov, and I. I. Oleynik, Phys. Rev. B 82, 064113 (2010). 23 A. B. Sivak, V. A. Romanov, and V. M. Chernov, Crystallogr. Rep. 55, 97 (2010). 104109-7 PHYSICAL REVIEW B 84, 104109 (2011) SERGEY V. STARIKOV et al. 24 N. Soneda, S. Ishino, and T. D. de la Rubia, Philos. Mag. Lett. 81, 649 (2001). 25 C. Becquart, C. Domain, A. Legris, and J. van Duysen, J. Nucl. Mater. 280, 73 (2000). 26 R. Pasianot, M. Alurralde, A. Almazouzi, and M. Victoria, Philos. Mag. A 82, 1671 (2002). 27 Y. N. Osetsky, D. J. Bacon, A. Serra, B. N. Singh, and S. Golubov, Philos. Mag. 83, 61 (2003). 28 J.-H. Shim, H.-J. Lee, and B. D. Wirth, J. Nucl. Mater. 351, 56 (2006). 29 A. F. Calder, D. J. Bacon, A. V. Barashev, and Y. N. Osetsky, Philos. Mag. Lett. 88, 43 (2008). 30 K. O. E. Henriksson, K. Nordlund, and J. Keinonen, Phys. Rev. B 76, 245428 (2007). 31 S. Han, L. A. Zepeda-Ruiz, G. J. Ackland, R. Car, and D. J. Srolovitz, Phys. Rev. B 66, 220101(R) (2002). 32 D. Nguyen-Manh, A. P. Horsfield, and S. L. Dudarev, Phys. Rev. B 73, 020101 (2006). 33 S. Rokkam, A. El-Azab, P. Millett, and D. Wolf, Modell. Simul. Mater. Sci. Eng. 17, 064002 (2009). 34 S. L. Dudarev, M. R. Gilbert, K. Arakawa, H. Mori, Z. Yao, M. L. Jenkins, and P. M. Derlet, Phys. Rev. B 81, 224107 (2010). 35 P. Ehrhart, P. Jung, H. Shultz, and H. Ullmaier, in Atomic Defects in Metals, edited by H. Ullmaier (Springer, Berlin, 1991), Vol. 25. 36 M. W. Finnis and J. E. Sinclair, Philos. Mag. A 50, 45 (1984). 37 P. M. Derlet, D. Nguyen-Manh, and S. L. Dudarev, Phys. Rev. B 76, 054107 (2007). 38 M. S. Daw and M. I. Baskes, Phys. Rev. B 29, 6443 (1984). 39 The presented Mo-Xe EAM potential can be sent by the authors via email upon request. 40 F. Ercolessi and J. B. Adams, Europhys. Lett. 26, 583 (1994). 41 P. Brommer and F. Gahler, Modell. Simul. Mater. Sci. Eng. 15, 295 (2007). 42 G. Kresse and J. Furthmuller, Phys. Rev. B 54, 11169 (1996). 43 J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, and C. Fiolhais, Phys. Rev. B 46, 6671 (1992). 44 J. P. Biersack and J. F. Ziegler, Nucl. Instrum. Methods Phys. Res. 194, 93 (1982). 45 K. Syassen and W. B. Holzapfel, Phys. Rev. B 18, 5827 (1978). 46 Y. K. Vohra and A. L. Ruoff, Phys. Rev. B 42, 8651 (1990). 47 J. W. Edwards, R. Speiser, and H. L. Johnston, J. Appl. Phys. 22, 424 (1951). 48 D. Alfe, M. J. Gillan, and G. D. Price, J. Chem. Phys. 116, 6170 (2002). 49 The script for calculation of the melting line (see at [http://lammps.sandia.gov/scripts.html]). 50 A. B. Belonoshko, S. I. Simak, A. E. Kochetov, B. Johansson, L. Burakovsky, and D. L. Preston, Phys. Rev. Lett. 92, 195701 (2004). 51 R. S. Hixson, D. A. Boness, J. W. Shaner, and J. A. Moriarty, Phys. Rev. Lett. 62, 637 (1989). 52 D. Errandonea, B. Schwager, R. Ditz, C. Gessmann, R. Boehler, and M. Ross, Phys. Rev. B 63, 132104 (2001). 53 C. Cazorla, M. J. Gillan, S. Taiol, and D. Alfe, J. Phys.: Conf. Series 121, 012009 (2008). 54 R. Boehler, D. Santamaria-Perez, D. Errandonea, and M. Mezouar, J. Phys.: Conf. Series 121, 022018 (2008). 55 S. J. Plimpton, J. Comput. Phys. 117, 1 (1995). 56 See Supplemental Material at http://link.aps.org/supplemental/ 10.1103/PhysRevB.84.104109 for illustration of the structure: the animation shows the structure of SIAC from 5 SIAs. 57 See Supplemental Material at http://link.aps.org/supplemental/ 10.1103/PhysRevB.84.104109 for illustration of the evolution of SIAs at temperature 1200 K. The 1D thermal diffusion of SIAC is shown. Total time of simulation is equal to 1 ns. Only the atoms with high potential energy are shown (atoms forming SIAs). 58 See Supplemental Material at http://link.aps.org/supplemental/ 10.1103/PhysRevB.84.104109 for movie: the animation shows the vacancy loop formation near surface at radiation damage. Only the disordered atoms (atoms forming defects) are shown. T = 1100 K. 59 F. Gao and D. J. Bacon, Philos. Mag. A 75, 1603 (1997). 60 S. Muller, M. L. Jenkins, C. Abromeit, and H. Wollenberger, Philos. Mag. A 75, 1625 (1997). 104109-8