Improvements to the Sink Strength Theory Used in Multi-Scale Rate Equation Simulations of Defects in Solids

Abstract

The application of mean-field rate theory equations have proven to be a versatile method in simulating defect dynamics and temporal changes in the microstructure of materials. The reliability and usefulness of the method, however, depends critically on the defect interaction parameters used. In this study, we show that the main interaction parameter, the sink strength, intrinsically depends on the detrapping, or the dissociation process itself. We present a theory on how to determine the appropriate sink strengths. The correct sink strength required for a detrapping defect, is considerably larger than the values commonly used, and thus should not be neglected.

1. Introduction

The physical and mechanical properties of materials are largely determined by their microstructure, and defect and impurity concentrations. To understand and control these changes during aging, ion irradiation or annealing, a simulation technique extending to long time and length scales is required. Among simulation techniques able to fulfil these challenging demands are the mean-field rate theory equations (RE) and the kinetic Monte Carlo (KMC) methods. KMC is a stochastic simulation method, where all the dynamic properties and reactions for all incorporating defects have to be known a priori. A single time step used in a KMC simulation, however, might be of the order of 10⁻¹² s for fast occurring processes [1], which will restrict the accessible time and the amount of defect concentrations studied with this method.

In the mean-field REs, all the relevant defect mobilities and reactions are collected to a set of non-linear partial differential equations that are temporally and spatially solved [2,3,4,5]. The RE methodology has been extensively used for simulating various dynamic processes in solid materials. These studies include trapping simulations of helium (He) and hydrogen (H) isotopes in vacancies [6,7], clustering of irradiation-induced vacancies and self-interstitials [3,8,9,10], He and H bubble formation [11,12], H diffusion and retention [13], precipitation, swelling, and segregation [14,15,16], fusion edge-localized plasma modes simulations [17], H isotope exchange [18,19,20], irradiation-induced amorphization [21], and simulation of thermal desorption spectrometry (TDS) profiles [22,23,24,25] to mention just a few.

The parameters for detrapping or dissociation processes can be determined by reaction rate theory [26] and transition state theory [27,28]. The trapping or annihilation processes are described by the so-called sink strength parameter [2,3,15]. Sink strengths have been determined analytically for various symmetric defects or traps including spherical traps, dislocation lines, and grain boundaries [3,14,29,30]. Finding the sink strength for arbitrarily shaped traps the conventional Monte Carlo (MC) method can be used [31,32,33].

The main advantage of the RE is that it is a very fast method for simulating processes taking place for extended time and length scales. The usefulness of the simulation results, however, depends critically on the defect interaction parameters, i.e., there is a critical requirement that the sink strengths applied, describe the physical processes correctly. In this paper, we show that the sink strength depends decisively on the detrapping or the dissociation process. We further derive the equations highlighting the requirements how the sink strengths and REs should be modified accordingly. Finally, we make test simulations to compare the results obtained using different sink strength theories. 2. Results 2.1. Theory

The presentation of the rate equation and sink strength theories are given in Refs. [2,3,14,15]. The governing idea is to relate the time evolution of any defect concentration to the trapping and detrapping processes. A simplified rate equation with only one kind of trap, excluding boundary diffusion and source terms for brevity, is given by

dcdt=−_ce∮^FindA+_cf∮^FoutdA=−Dkc+E,

where_ceand_cfare the concentrations of the empty and filled traps,^Finand^Foutare the trapping and detrapping fluxes of defect c to and from a trap, respectively. Parameters_ce∮^FindAand_cf∮^FoutdAdenote the trapped and detrapped defect concentrations per second, where A is the trap surface area. These integrals are impractical to use and are substituted in the RE method by terms−Dkcand+E , respectively, as presented in Equation (1). D is the defect diffusion coefficient, k the sink strength, and E is the thermal emission (detrapping) of defects c from the filled traps_cf. According to rate theory, E can be represented as

E=_cfνexp−_EtkBT,

whereνis the detrapping attempt frequency,_Et(≈_Eb+_Em) is the trapping or dissociation energy,_Ebis the binding energy and_Emthe defect migration barrier,_kBis the Boltzmann constant, and T the absolute temperature.

Usually, the sink strengths, k, are determined analytically or by applying the MC method. The detrapping attempt frequency and the trapping energy are obtained with atomistic simulations or with experiments. As an example of analytical sink strength, and to be used later in this work, the widely applied sink strength for spherical traps is given by Brailsford and Bullough [29]

k=4π_{Rt ce}(1+_Rtk),

where_Rtis the trapping radius and_ce the trap concentration. Most of the sink strengths used in the literature are defined for randomly distributed defects in the material which result, for instance, from ion irradiation of the material. The same sink strengths are usually used irrespective of how the defects enter the material. We will, however, now show that the sink strengths influencing detrapped defects, are greatly enhanced compared to sink strengths of defects entering the material through irradiation. This is done by comparing sink strengths for defects having the initial position either randomly distributed in the cell, or in close proximity of a trap boundary. The fast MC method [33] was used in calculating the sink strengths k for different defect positions with respect to the trap. The results are shown in Figure 1 as a function of the relative amount of traps, i.e., the trap volume fraction.

In the MC calculations, the traps simulated are spherical with a trapping radius of_Rt = 1.0 nm and the trap concentration is varied between 6 × 10⁻⁶ − 0.4 nm⁻³. The defect jump length is set to 0.1 nm. The number of defects for each trap concentration is 10 × 10⁶ giving a statistical error of about 0.03% for the sink strengths [33]. The results show that defects with an initial position close to a trap have a large probability to be trapped in the adjacent trap. Thus, their sink strengths are much larger than the sink strengths for defects with the random position further away from any trap. In addition, in Figure 1 are shown the analytical sink strength values as obtained from Ref. [33] and compared against the values obtained using Equation (3). The analytical results show minor deviation features only at higher trap volume fractions, which will be discussed later in this work.

To formulate the aforementioned position-sensitive effect of sink strengths, we present the modified sink strength, which incorporates the position of the defect with respect to the trap. We introduce a dimensionless enhancement factor_ϵkas_ϵk=_kclose/_krand , which is defined as the sink strength for defects located close to the trap divided by the sink strength for defects with random initial position. Figure 2 shows this enhancement factor as a function of trap volume fraction for three different trapping radii: 0.5, 1.0, and 2.0 nm. The results for trapping radius 1.0 nm are calculated from the MC sink strength results presented in Figure 1.

It can be seen that the enhancement factor increases with the trap size, which reflects the fact that the trapping probability for the defect close to the trap increases as the trap size increases. The enhancement factor approaches a maximum value as the trap volume fraction approaches zero. As the trap volume fraction increases, the trap concentration becomes higher, and the probability for the defect to be trapped to another trap than the nearest one increases, leading to a smaller enhancement factor. For very high trap volume fractions the_ϵkapproaches unity, when the sink strengths_kcloseand_krand approach equal values. The key point is that the sink strengths for defects vary depending on their lattice position. These positions may be randomly distributed as a result of irradiation, or close to an empty trap due to detrapping. With this in mind, we revisit the RE theory as presented by Equation (1). At equilibrium (dc/dt=0 ) the trapping and detrapping rates are in balance and Equation (1) yields

Dkc=E=_ce∮^FindA=_{ϵk cf}∮^FoutdA,

which takes into account the sink strength enhancement_ϵk. From this we get

_cf∮^FoutdA=E/_ϵk,

which together with the sink strength from Equation (1) yields

k=1Dc_ce∮^FindA−E/_ϵk+E=_ceDc∮^FindA+E(1−1/_ϵk)Dc=_krand+E(1−1/_ϵk)Dc.

The first term in Equation (6),_krand , is the usual random position sink strength [3,15]. The second term contains the factor_ϵk, which enhances the sink strength when detrapping is activated. With no detrapping, or if the enhancement factor is one, the sink strength reduces to the one conventionally used. The full rate equation, including diffusion and the source of defects (S) originating from the irradiation process, now reads

dcdt=D^∇2c−Dkc+E+S,

where the k and E are given as in Equations (2) and (6), respectively. An alternative formulation that gives the right equilibrium condition is obtained when Equation (6) is inserted in Equation (7)

dcdt=D^∇2c−D_krandc+E/_ϵk+S.

2.2. Rate Equation Simulations

To highlight the effect of using different sink strength theories, we perform the following RE simulation for detrapping, trapping, and diffusion processes during annealing. The simulation conditions have been chosen such, that the corresponding KMC simulation is possible to be carried out allowing comparison. In these simulations, we have a 200 nm thick material layer with a Gaussian-shaped trap profile located at a mean depth of 50 nm from the surface and with a standard deviation (SD) 10 nm. The maximum trap concentration value at the mean depth is set to 1 × 10¹⁸ cm⁻³. Initially, the traps are filled with one defect per trap, and the trapping radius R_t is set to 1.0 nm. The defect detrapping attempt frequencyν is 5 × 10¹² Hz and the binding energy E_b = 0.8 eV, as used in Equation (2). The temperature at the beginning of the simulation is set to 300 K, and it increases linearly with a rate of 50 K/s up to 800 K during the 10 s simulation. As the temperature is rising in the simulation, detrapping, and diffusion processes of defects start to take place. Some defects are retrapped, but some reach the surface and leave the simulation layer. This flux of defects to the surface is monitored and the corresponding results obtained with RE and KMC methods are compared. This simulation represents the popular TDS method frequently used in the experiments and the corresponding computational analyses. Processes taking place on the surface have been omitted in the present simulations. Thus, the defect flux to surface is taken directly as the flux of defects crossing the simulation layer boundary. The defect diffusion parameters used are: jump length is 0.1 nm, jump frequency is 5 × 10¹² Hz, and the migration barrier_Em= 0.25 eV. The detrapping distance for the defects from the traps is chosen as 0.05 nm. In this simulation, the RE is represented in 1D and reads

dcdt=D^d2cd^z2−Dkc+E,

with z being the depth coordinate. The different sink strengths k and detrapping terms E used are given in Table 1. The difference in using Equation (3) or the sink strengths from Ref. [33] for_krand in the table is very small, as seen in Figure 1, due to the low trap volume fraction used in the simulation.

Figure 3 illustrates the resulted front surface fluxes of defects as obtained with the RE and KMC simulations. It is apparent from the results that the simulated TDS peak (i.e., the calculated defect flux through the front surface) occurs too fast with the usual RE formalism. In view of this study, this is to be expected, because the usually used random sink strength does not take into account the enhanced retrapping due to the defects being detrapped close to the empty traps. The results obtained with the new theory Equation (7) and the alternative formulation Equation (8) agree nicely with the KMC simulation. Moreover, the widely used formulation for sink strength and detrapping, (Table 1) [29], results in a delayed TDS release peak.

Equation (7) with the actual sink strength (given in Equation (6)) should be used if the correct trapped and detrapped defect rates (defect/s) are wanted from the simulation. In this case, the KMC and RE give the same results for the individual trapped and detrapped number of defects per second. Looking at the change in the defect concentration (c) as a function of time, however, both Equations (7) and (8) give the same result. The term (E/_ϵk ) in Equation (8) is actually the detrapped defects minus the retrapped defects due to the enhanced trapping. The reason for using Equation (8), is that it is simpler and that the time steps in the simulation can be chosen bigger.

This example underlines the importance of the presented theory, i.e., not using enhanced sink strengths, which take into account the defect lattice position, can lead to totally different results. In many studies, the different trapping parameters are fitted so that the RE simulation agrees with the experiment. With the wrong sink strength theory, these fitted parameters are questionable and may have no scientific value.

Further, in Figure 4 are shown snapshots of the depth profiles of traps and defects at different times during the performed RE and KMC simulations. The chosen times are the beginning, t = 0 s (Figure 4a), at t = 3 s (450 K) where detrapping has started (Figure 4b), at t = 4 s (500 K) where the defect flux through the front surface is at its maximum (Figure 4c), and at t = 5 s (550 K) where almost all traps have been emptied and the defect flux through the front surface approaches zero (Figure 4d).

The profiles are named as follows: traps with defects (Trap-Def), empty traps (Trap), and mobile free defects (Def). From Figure 4 we can see the RE and KMC simulations agreeing perfectly throughout the whole comparable concentration region. The computing time, however, limits the number of traps and defects that can be studied with the KMC method. This can be seen as the lower concentration limit of about 1 × 10¹⁴ defects/cm³ shown in Figure 4d), whereas the corresponding maximum concentration of the diffusing defect profile, obtained with the RE simulation, is well below 1 × 10¹⁴ cm⁻³. Consequently, it is not possible to deduce the temporal defect depth profiles from the KMC simulation. The perfect agreement of the TDS peaks obtained with the RE and KMC methods (Figure 3), however, indicates that the mobile defect profiles will continue to match at lower concentrations, if the number of defects in the KMC simulation would be sufficient.

It must be mentioned that even though the simulation set-up is chosen favourable for the KMC technique, the RE simulation is about five orders of magnitude faster than the KMC simulation (about 10 s compared to 11 days, respectively). The fact that the number of free defects in the KMC simulation is too small to be used in obtaining any useful defect depth profiles is a serious shortcoming. Often it is very crucial to examine the depth profiles in order to understand the different aspects of the experimental findings, and in obtaining other wanted physical quantities from the simulations. 3. Discussion The results of this study show that the sink strengths depend critically on the position-related detrapping or dissociation processes. Therefore, the numerous simulations including detrapping, done for about 50 years are probably not entirely accurate. The new sink strength formulation needs the new enhancement factor parameter to be determined, which is obtained using, for instance, the MC method. Including the detrapping effect in the sink strength interaction parameters will increase the usefulness of the RE method. The RE simulations should now give more reliable defect and trap parameters, and the long length and time scale microstructure simulations should be more accurate. 4. Materials and Methods The authors declare that the data supporting this study and the MC and RE codes are available from the corresponding author upon request.

Figure 1. The sink strengths for spherical traps calculated by Monte Carlo (MC) and analytical methods. The position for the defect in the simulation cell is either random or at a close distance (0.05 nm) from a trap boundary. The sink strength for defects with close to trap initial position is greatly enhanced compared to the ones with the random initial position.

Figure 2. The enhancement factorϵkfor sink strengths as a function of trap volume fraction for three different trapping radii. At high trap concentrations, the probability of a defect to be trapped to other than the nearest trap increases, which is seen as a decrease in theϵkvalues.

Figure 3. The front surface defect flux comparison between the kinetic Monte Carlo (KMC) and rate theory equations (RE) simulations using different sink strength theories (see text and Table 1).

Figure 4. The depth profiles for the KMC and RE simulations at four different times and temperatures during the TDS simulation. At the beginning, at time 0 s, (a) the temperature is 300 K, and there is only the filled trap profile (Trap-Def). The temperature increases linearly with time (50 K/s) and at 3 s, (b) some defects (Def) have detrapped, leaving empty traps (Trap). At 4 s, (c) there are more empty than filled traps, and at 5 s, (d) only a few filled traps are left. Note that the number of defects (Def) for the KMC simulation are too few to determine any meaningful defect depth profiles.

	Sink Strength (k)	Detrapping (E)
(1) Usual	krand	cfνexp(−Et/kBT)
(2) New	krand+E(1−1/ϵk)/(Dc)	cfνexp(−Et/kBT)
(3) New (Alt) a	krand	cfνexp(−Et/kBT)/ϵk
(4) B&B b	krand	4πRt cf(1+Rtk)Dexp(−Eb/kBT)

a Alternative formulation Equation (8). b Brailsford and Bullough [29].

Author Contributions

Conceptualization, T.A.; methodology, T.A.; software, T.A.; validation, T.A.; formal analysis, T.A.; investigation, T.A.; writing-original draft preparation, T.A. and K.H.; writing-review and editing, T.A. and K.H.; visualization, T.A. All authors have read and agreed to the published version of the manuscript.

Funding

This work has been carried out within the framework of the EUROfusion Consortium, and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

Acknowledgments

Computational resources provided by the Finnish Grid and Cloud Infrastructure (persistent identifier urn:nbn:fi:research-infras-2016072533) is gratefully acknowledged.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

1. Derlet, P.M.; Nguyen-Manh, D.; Dudarev, S.L. Multiscale modeling of crowdion and vacancy defects in body-centered-cubic transition metals. Phys. Rev. B 2007, 76, 054107.

2. McNabb, A.; Foster, P.K. A new analysis of the diffusion of hydrogen in iron and ferritic steels. Trans. Metall. Soc. AIME 1963, 227, 618.

3. Wiedersich, H. On the theory of void formation during irradiation. Radiat. Eff. 1972, 12, 111-125.

4. Freeman, G. Kinetics of Nonhomogeneous Processes, LK Mansur: Mechanisms and Kinetics of Radiation Effects in Metals and Alloys; John Wiley and Sons: New York, NY, USA, 1987.

5. Ahlgren, T.; Heinola, K.; Vörtler, K.; Keinonen, J. Simulation of irradiation induced deuterium trapping in tungsten. J. Nucl. Mater. 2012, 427, 152-161.

6. Baskes, M.I.; Wilson, W.D. Kinetics of helium self-trapping in metals. Phys. Rev. B 1983, 27, 2210-2217.

7. Myers, S.M.; Nordlander, P.; Besenbacher, F.; Nørskov, J.K. Theoretical examination of the trapping of ion-implanted hydrogen in metals. Phys. Rev. B 1986, 33, 854-863.

8. Katoh, Y.; Muroga, T.; Kohyama, A.; Stoller, R.E.; Namba, C. Rate theory modeling of defect evolution under cascade damage conditions: The influence of vacancy-type cascade remnants on defect evolution. J. Nucl. Mater. 1996, 233-237, 1022-1028.

9. Rottler, J.; Srolovitz, D.J.; Car, R. Point defect dynamics in bcc metals. Phys. Rev. B 2005, 71, 064109.

10. Ortiz, C.J.; Caturla, M.J. Simulation of defect evolution in irradiated materials: Role of intracascade clustering and correlated recombination. Phys. Rev. B 2007, 75, 184101.

11. Wilson, W.D.; Baskes, M.I.; Bisson, C.L. Atomistics of helium bubble formation in a face-centered-cubic metal. Phys. Rev. B 1976, 13, 2470-2478.

12. Wang, Z.; Sang, C.; Liu, S.; Chang, M.; Sun, J.; Wang, D. Modeling of fuel retention in the pre-damaged tungsten with MeV W ions after exposure to D plasma. Nucl. Mater. Energy 2017, 13, 1-7.

13. Coenen, J.; Berger, M.; Demkowicz, M.; Matveev, D.; Manhard, A.; Neu, R.; Riesch, J.; Unterberg, B.; Wirtz, M.; Linsmeier, C. Plasma-wall interaction of advanced materials. Nucl. Mater. Energy 2017, 12, 307-312.

14. Wert, C.; Zener, C. Interference of Growing Spherical Precipitate Particles. J. Appl. Phys. 1949, 21, 5.

15. Brailsford, A.; Bullough, R. The rate theory of swelling due to void growth in irradiated metals. J. Nucl. Mater. 1972, 44, 121-135.

16. Barnard, L.; Tucker, J.; Choudhury, S.; Allen, T.; Morgan, D. Modeling radiation induced segregation in Ni-Cr model alloys from first principles. J. Nucl. Mater. 2012, 425, 8-15.

17. Heinola, K.; Ahlgren, T.; Brezinsek, S.; Vuoriheimo, T.; Wiesen, S. Modelling of the effect of ELMs on fuel retention at the bulk W divertor of JET. Nucl. Mater. Energy 2019, 19, 397-402.

18. Schmid, K.; von Toussaint, U.; Schwarz-Selinger, T. Transport of hydrogen in metals with occupancy dependent trap energies. J. Appl. Phys. 2014, 116, 134901.

19. Hodille, E.A.; Ferro, Y.; Fernandez, N.; Becquart, C.S.; Angot, T.; Layet, J.M.; Bisson, R.; Grisolia, C. Study of hydrogen isotopes behavior in tungsten by a multi trapping macroscopic rate equation model. Phys. Scr. 2016, 2016, 014011.

20. Markelj, S.; Založnik, A.; Schwarz-Selinger, T.; Ogorodnikova, O.; Vavpetič, P.; Pelicon, P.; Čadež, I. In situ NRA study of hydrogen isotope exchange in self-ion damaged tungsten exposed to neutral atoms. J. Nucl. Mater. 2016, 469, 133-144.

21. Swaminathan, N.; Morgan, D.; Szlufarska, I. Ab initio based rate theory model of radiation induced amorphization in β-SiC. J. Nucl. Mater. 2011, 414, 431-439.

22. Pisarev, A.A.; Voskresensky, I.D.; Porfirev, S.I. Computer modeling of ion implanted deuterium release from tungsten. J. Nucl. Mater. 2003, 313-316, 604.

23. Poon, M.; Haasz, A.; Davis, J. Modelling deuterium release during thermal desorption of D+-irradiated tungsten. J. Nucl. Mater. 2008, 374, 390-402.

24. Gasparyan, Y.; Ogorodnikova, O.; Efimov, V.; Mednikov, A.; Marenkov, E.; Pisarev, A.; Markelj, S.; Čadež, I. Thermal desorption from self-damaged tungsten exposed to deuterium atoms. J. Nucl. Mater. 2015, 463, 1013-1016.

25. Grigorev, P.; Matveev, D.; Bakaeva, A.; Terentyev, D.; Zhurkin, E.E.; Oost, G.V.; Noterdaeme, J.M. Modelling deuterium release from tungsten after high flux high temperature deuterium plasma exposure. J. Nucl. Mater. 2016, 481, 181-189.

26. Arrhenius, S. Über die Reaktionsgeschwindigkeit bei der Inversion von Rohrzucker durch Säuren. Z. Phys. Chem. (Leipzig) 1889, 4, 226.

27. Wigner, E. On the penetration of potential energy barriers in chemical reactions. Z. Phys. Chem. Abt. B 1932, 19, 203.

28. Eyring, H. The activated complex in chemical reactions. J. Chem. Phys. 1935, 3, 107.

29. Brailsford, A.D.; Bullough, R. Theory of sink strengths. Philos. Trans. R. Soc. Lond. 1981, 302, 87.

30. Rouchette, H.; Thuinet, L.; Legris, A.; Ambard, A.; Domain, C. Quantitative phase field model for dislocation sink strength calculations. Comput. Mater. Sci. 2014, 88, 50-60.

31. Malerba, L.; Becquart, C.S.; Domain, C. Object kinetic Monte Carlo study of sink strengths. J. Nucl. Mater. 2007, 360, 159.

32. Jansson, V.; Malerba, L.; Backer, A.D.; Becquart, C.S.; Domain, C. Sink strength calculations of dislocations and loops using OKMC. J. Nucl. Mater. 2013, 442, 218-226.

33. Ahlgren, T.; Bukonte, L. Sink strength simulations using the Monte Carlo method: Applied to spherical traps. J. Nucl. Mater. 2017, 496, 66.