The role of Frenkel defect diffusion in dynamic annealing in ion-irradiated Si

The formation of stable radiation damage in crystalline solids often proceeds via complex dynamic annealing processes, involving migration and interaction of ballistically-generated point defects. The dominant dynamic annealing processes, however, remain unknown even for crystalline Si. Here, we use a pulsed ion beam method to study defect dynamics in Si bombarded in the temperature range from −20 to 140 °C with 500 keV Ar ions. Results reveal a defect relaxation time constant of ~10–0.2 ms, which decreases monotonically with increasing temperature. The dynamic annealing rate shows an Arrhenius dependence with two well-defined activation energies of 73 ± 5 meV and 420 ± 10 meV, below and above 60 °C, respectively. Rate theory modeling, bench-marked against this data, suggests a crucial role of both vacancy and interstitial diffusion, with the dynamic annealing rate limited by the migration and interaction of vacancies.

The formation of stable radiation damage in crystalline materials often proceeds via so-called dynamic annealing (DA) processes, involving migration, recombination, and clustering of mobile point defects during irradiation. These DA processes are complex and remain poorly understood even for crystalline Si, which is the most extensively studied and arguably best understood material 1 . Since DA is believed to be thermally activated, determining activation energies (E a s) of the dominant DA processes is one of the first logical steps toward understanding radiation damage dynamics.
Several attempts to understand DA in Si by measuring associated E a s have been reported, albeit with limited success as they have revealed a very wide range of E a s of ~0.2-1.7 eV [2][3][4][5] . For example, Linnros and Holmen 2 have extracted an E a of 1.2 eV from the temperature (T) dependence of the ion dose rate at which the crystalline/amorphous interface is stationary for bombardment with 1.5 MeV Ne, Ar, or Xe ions in the T range of ~100-300 °C. In contrast, Schultz et al. 3 have found an E a of 0.9 eV by plotting the T-dependence of the dose rate required to reach amorphization in the crystal bulk at a fixed dose for 1 MeV Si ion irradiation in a narrow T range of ~40-80 °C. Goldberg et al. 4 have expanded the approach of Schultz et al. 3 and found E a s in a wide range of ~0.7-1.7 eV for 80 keV ions with masses ranging from 12 C to 132 Xe and T in the range of ~10-300 °C. They 4 have found that the E a value increases close-to-linearly with either ion mass or T (since, with increasing ion mass, the onset of bulk amorphization occurs at higher Ts for any given dose and dose rate). Finally, Kinomura et al. 5 , following the work of Linnros et al. 6 , have systematically studied T dependencies of the rate of ion-beam-induced epitaxial crystallization and found E a values of ~0.3-0.4 eV in the T range of ~250-400 °C for irradiation with 3 MeV C, Si, Ge, or Au ions and an E a of ~0.2 eV for 3 MeV C ion bombardment in the T range of ~150-280 °C. Such large inconsistency in the E a s reported (~0.2-1.7 eV) highlights the complexity and currently limited understanding of the fundamental physics governing DA even for Si.
Here, we use a novel pulsed ion beam technique 7-10 to measure the T dependence of the DA time constant (τ) in Si bombarded with 500 keV Ar ions in a regime of relatively high ion doses when damage accumulation is dominated by inter-cascade DA processes (i.e., by the interaction of mobile defects generated in different collision cascades). Our results reveal two well-defined regions in the Arrhenius plot of the DA rate, with very different E a s of 73 and 420 meV, below and above ~60 °C, respectively. A comparison of these E a s with the literature values [11][12][13][14][15] and our rate theory modeling results suggest that inter-cascade communication in Ar-ion-bombarded Si is carried out primarily by migrating interstitials and vacancies, with vacancy migration and interaction being the rate limiting processes.
Scientific RepoRts | 7:39754 | DOI: 10.1038/srep39754 Experimental The 4 MV ion accelerator (National Electrostatics Corporation, model 4UH) at Lawrence Livermore National Laboratory was used for both ion irradiation and ion beam analysis. Float-zone grown (100) Si single crystals (with a resistivity of ~5 Ω cm) were bombarded with 500 keV 40 Ar + ions at 7° off the [100] direction in the T range from − 20 to 140 °C. To improve thermal contact, the samples were attached to the Cu sample holder with conductive Ag paste. All irradiations were performed in a broad beam mode 7 . In each irradiation run, the total dose was split into a train of equal square pulses each with an instantaneous dose rate F on ≈ 1.9 × 10 13 cm −2 s −1 and duration t on = 1 ms, corresponding to ~4.6 × 10 −5 displacements per atom per pulse. The depth profile of ballistically-generated vacancies was calculated with the TRIM code (version SRIM-2013.00) 16 with an atomic concentration of Si of 4.98 × 10 22 atoms cm −3 and a threshold energy for atomic displacements of 13 eV. The pulsing parameters F on and t on were chosen based on previous pulsed beam measurements of Si at room T 7,8 in order to maximize the DA efficiency and to limit the inter-cascade defect interaction within each pulse. The adjacent pulses were separated by time t off , which was varied between 0.2 and 50 ms. The inset in Fig. 1(a) shows a schematic of the time dependence of the instantaneous dose rate and defines pulsing parameters t on , t off , and F on . A more detailed description of the experimental arrangement can be found elsewhere 7,8 .
The dependence of lattice damage on t off was studied ex-situ at room T by ion channeling. Depth profiles of lattice disorder were measured with 2 MeV 4 He + ions incident along the [100] direction and backscattered into a detector at 164° relative to the incident beam direction. Raw channeling spectra were analyzed with one of the conventional algorithms 17 for extracting depth profiles of relative disorder. Values of average relative bulk disorder (n) were obtained by averaging depth profiles of relative disorder over 20 channels (~38 nm) centered on the bulk damage peak maximum. Error bars of n are standard deviations. Total ion doses at different Ts were different and chosen such that, for continuous beam irradiation (i.e., t off = 0), n was in the range of 0.6-0.9 (with n = 1 corresponding to full amorphization). Results and Discussion Figure 1 shows representative depth profiles of relative disorder for bombardment of Si with continuous (t off = 0 ms) and pulsed (t off = 3 and 30 ms) beams at Ts of 40 and 80 °C. It is seen that, for both Ts, these depth profiles are bimodal, with the first small peak at the sample surface and the second major peak in the crystal bulk. The bulk peak is centered on ~500 nm, which corresponds to the maximum of the nuclear energy loss profile for 500 keV Ar ions 16 . It is seen from Fig. 1 that the average bulk disorder (n) decreases with increasing t off for both Ts. This is better illustrated in Fig. 2, which summarizes all the n(t off ) dependencies for Ts from − 20 to 140 °C. It is seen from Fig. 2 that, for all the Ts studied, n monotonically decreases with increasing t off . This effect is due to the interaction of mobile defects generated in different pulses and, hence, in different cascades (i.e., inter-cascade defect interaction).
Decay dependencies, as revealed by Fig. 2, are commonly described by either first or second order decay equations in order to estimate decay time constants and, hence, kinetic rates. In a first kinetic order process, the rate is directly proportional to the concentration of interacting species, while in the second order process, the rate is proportional to the square of the concentration. Examples of the first order processes of defect interaction include trapping of interstitials or vacancies at sinks, whereas direct vacancy-interstitial annihilation and the formation of di-vacancies and di-interstitials are examples of the second order processes.
The n(t off ) dependencies from Fig. 2 were fitted via the Marquardt-Levenberg algorithm 18 with first (n(t off ) = n ∞ + (n(0) − n ∞ ) exp(− t off /τ 1 )) and second order   Here, τ 1,2 is the characteristic decay time constant, and n ∞ is relative disorder for τ  t off 1,2 . Since a decrease in n with increasing t off is a result of inter-cascade defect interaction, the time constant τ reflects inter-cascade defect interaction processes. Below 60 °C, the data is fitted best with the second order decay equation, while the first order decay gives a better fit for all the n(t off ) dependencies above 60 °C. Such an evaluation of the kinetic order of n(t off ) dependencies was done by comparing R-squared values of fits with the first and second order decay equations. In all the cases of different Ts, however, R-squared values were > 0.96.
Temperature dependencies of τ 1 and τ 2 are plotted in Fig. 3, revealing a monotonic decrease with increasing T and a kink (i.e., a change in the first derivative) at 60 °C in both curves. As expected from the form of first and second order decay equations, τ 1 > τ 2 . Also plotted in Fig. 3 is the T dependence of the DA efficiency ξ 1,2 , which we define as before refs [7,8]: ξ = (n(0) − n ∞ )/n(0). Figure 3 shows that, within experimental errors, ξ increases with T, reflecting a corresponding decrease in n ∞ , which is also clearly seen in n(t off ) dependencies of Fig. 2. Above ~60 °C, ξ 1 saturates at ~90%. Note that the apparent saturation of ξ 2 at 100% for higher Ts is an artifact of inferior fitting with the second order decay equation above 60 °C.
Values of E a s obtained by fitting with the first and second order decay equations are comparable. Indeed, fitting data from Fig. 4 with the first order decay equation for T < 60 °C and with the second order decay equation for  T > 60 °C gives E a s and As of E a1 = 110 ± 10 meV and A 1 = (15.3 ± 5.4) × 10 3 Hz for T < 60 °C and E a2 = 440 ± 20 meV and A 2 = (4.05 ± 2.9) × 10 9 Hz for T > 60 °C.
Our experimental data unambiguously reveals the existence of two distinct dominant DA process at Ts below and above 60 °C, evidenced by (i) the presence of two well-defined Arrhenius regions with vastly different E a s and As (Fig. 4), (ii) a saturation of ξ for T > 60 °C (Fig. 3), and (iii) the change from the second order to the first order kinetic behavior. The switch from the second kinetic order to the first order behavior at 60 °C points to a change in the dynamics of defect interaction. Our Raman scattering measurements (with 633 nm laser light) of samples irradiated at different Ts, however, have not revealed any evidence of the change in the damage state at 60 °C. In addition, our transmission electron microscopy analysis has not revealed any differences in the type of damage in Si irradiated with pulsed and continuous beams to the same disorder level.
The two E a s of ~73 and ~420 meV (Fig. 4) are in contrast to much larger E a values of ~0.7-1.7 eV reported in previous dose rate studies [2][3][4] . In most of these prior attempts to measure the E a of DA 2-4 , the dose rate was used as the rate of the kinetic process in the Arrhenius relationship, with the implicit assumption that the dose rate is proportional to the rate of the dominant DA process. The large difference between the E a values measured in the present work and those in previous dose rate studies [2][3][4] could be attributed to the fact that, in the dose rate approach, the E a is effectively extracted from the ξ(T) dependence. As discussed in detail recently 8 , for our choice of F on and t on , ξ is the magnitude of the dose rate effect; i.e., the difference between n for continuous beam irradiation with dose rates of F = F on and F → 0. Hence, ξ reflects the fraction of ballistically-generated Frenkel defects that participate in DA processes for any given F on rather than the rate of defect interaction. In other words, while τ reflects the DA rate (i.e., dynamics), ξ describes the DA "magnitude". This is further supported by Fig. 3, revealing qualitatively different ξ(T) and τ(T) dependencies, indicating that parameters ξ and τ provide complementary information.
With such a large spread in the E a values reported previously [2][3][4][5] , the dominant DA processes in Si have remained elusive, with suggestions that the mobile defects in the high-dose regime dominated by inter-cascade DA processes behave differently from the migration of vacancies and interstitials in the low-dose regime of intra-cascade DA effects, commonly monitored by electron paramagnetic resonance (EPR) 11 , positron annihilation spectroscopy (PAS) 12 , and deep level transient spectroscopy (DLTS) 15 . In contrast, our results strongly suggest that, even at relatively high doses and dose rates typical for technologically relevant radiation environments, point defect migration still plays a key role in inter-cascade DA. We first note that an E a of ~73 meV for T < 60 °C agrees with E a s previously associated with the migration of isolated interstitials in Si revealed in DLTS (ref. 15) and PAS (ref. 12) measurements. Similarly, an E a of 420 meV for T > 60 °C agrees with E a s assigned to the migration of neutral vacancies (~0.40-0.46 eV) reported in EPR 11 , PAS 12 , and Raman spectroscopy 13 measurements. It is also consistent with several theoretical predictions [19][20][21][22][23] . We note that charge states of mobile defects are important for low dose irradiation conditions when lattice dopants determine the position of the Fermi level and the defect charge state. In contrast, the present study focuses on the regime of relatively high doses (and, hence, high damage levels) when we are likely dealing with only neutral vacancies and interstitials. Indeed, in the present study, the concentration of radiation-generated (stable) lattice defects largely exceeds the initial dopant concentration, and the material is in a semi-insulating state 8 . This conclusion is further supported by a recent pulsed beam study which showed τ to be independent of the doping level for 500 keV Ar ion irradiation of Si at RT 8 .
The dominant role of defect migration in inter-cascade DA is also consistent with a relatively large average separation between the centers of adjacent cascades in every pulse (~72 nm or ~300 atomic spacings). In such cases, inter-cascade defect interaction could not proceed without significant defect migration. For a three-dimensional random walk with a defect moving only a single atomic spacing during each jump, ~200,000 jumps would be necessary for the defect to travel the distance between the centers of two adjacent collision cascades 24 . Hence, while the E a of defect migration may not be the largest energy barrier a defect on the path to recombination or trapping must overcome, it is a barrier which must be overcome numerous times.
To get insight into the atomistics of defect interaction and to better correlate the E a s measured here with energetic barriers of specific defect migration or interaction processes, we have implemented the rate theory modeling as in refs 25-28, bench-marked against our experimental data. A successful description of both n(t off ) and τ(T) dependencies has been obtained with a model considering only the following four processes: (i) ballistic Frenkel pair generation (Is and Vs), calculated with the TRIM code 16 ; (ii) the annihilation of Is and Vs at unsaturating sinks (such as point defect clusters), (iii) V + V → V 2 , and (iv) I + V 2 → V. Here, I, V, and V 2 refer to interstitials, vacancies, and divacancies, respectively. The V + I → ∅ reaction is omitted since it does not affect the balance between V and I, which is critical for stable damage accumulation within this model. We have also found that adding the V + I → ∅ reaction to the model has no effect on the resultant E a values. The equations for T-dependent interaction parameters were taken from ref. 28. Migration energies of Vs and Is were set to 400 and 100 meV, respectively, while all other energetic barriers were set to zero. All capture radii were set to 5λ, where λ is the atomic spacing in Si. The total dose and the dose rate were set to 10 12 cm −2 and 10 13 cm −2 s −1 , respectively, with a t on of 1 ms. The attempt frequency and the sink concentration were set to 10 −11 s −1 and 8.5 × 10 16 cm 3 , respectively, and the bulk disorder (n) was represented by the V 2 concentration, as in previous studies [25][26][27][28] .
Results of such rate theory modeling are shown in Fig. 4, along with linear fits (dashed lines) to determine E a s. It is remarkable that, such a relatively simple model is capable of reproducing a rather complex experimental data set: n(t off ) dependencies and the Arrhenius behavior of the DA rate. Within this model, an E a of ~400 meV in the Arrhenius plot of Fig. 4 for T > 60 °C indeed corresponds to the V migration energy. However, an E a of ~100 meV for T < 60 °C in Fig. 4 actually does not correspond to the I migration energy. Instead, it reflects a competition between the two channels for V annihilation, one by trapping at sinks and the other via the formation of V 2 . While V 2 formation dominates at lower Ts, V annihilation at sinks becomes the rate limiting process at T > 60 °C. The much lower E a for I migration results in an I migration rate that is orders of magnitude faster than that of Vs. As a Scientific RepoRts | 7:39754 | DOI: 10.1038/srep39754 result, τ is controlled predominantly by the V migration rate in the entire T range studied. Hence, modeling also helps us understand the origin of the critical T of 60 °C in Fig. 4.
However, this relatively simple model has limitations. It takes into account only a very limited subset of possible defect interactions and ignores possible contributions from interstitial clusters and larger vacancy clusters. This simple model cannot quantitatively describe the full range of the damage buildup up to lattice amorphization as it does not take into account the non-linear processes leading to a super-linear (sigmoidal) damage accumulation at elevated Ts. This model also does not include nonlinear cascade density effects leading to the formation of thermal and/or displacement spikes 29 . In addition, this model does not predict the switch from the 2nd to the 1st order kinetic behavior at 60 °C. More work is currently needed for the development of more sophisticated and physically realistic models that can describe the entire range of experimental observations, including the damage buildup and defect interaction dynamics under different irradiation conditions. Future experimental work is also needed to study non-linear spike effects on the temperature dependence of defect interaction dynamics in Si.

Conclusion
In conclusion, we have used the pulsed beam technique to measure the T-dependence of the effective time constant of DA in Si bombarded in the T range from − 20 to 140 °C with 500 keV Ar ions. Results have revealed two well-defined Arrhenius regimes described by activation energies of ~73 and ~420 meV, below and above 60 °C, respectively. A comparison of experimental data with rate theory modeling results suggests that inter-cascade DA in Si proceeds via interstitial and vacancy migration, and the DA rate is limited by processes of vacancy annihilation at sinks and divacancy formation. Our results could have far reaching implications for future studies of radiation defect dynamics as the pulsed-beam method offers a novel way of measuring DA rates and activation energies that could be directly compared with results of predictive modeling. Importantly, this method can be used to measure defect interaction rates under technologically relevant irradiation conditions of pronounced inter-cascade defect interaction in any material, opening the door to studying intricate defect interaction phenomena in materials with the composition across the periodic table.