Kinetic Monte Carlo models and simulations have been employed to study cascade ageing and defect accumulation at low doses using the input from atomistic simulations of cascades and defect migration properties (Barashev, Bacon et al. 1999; Gao, Bacon et al. 1999; Barashev, Bacon et al. 2000; Heinisch, Singh et al. 2000; Heinisch and Singh 2003; Domain, Becquart et al. 2004; Becquart, Domain et al. 2005; Caturla, Soneda et al. 2006). These have mostly focused on intrinsic defect (interstitial and vacancy) migration, clustering and annihilation under irradiation conditions. In this example, we focus on extrinsic gas atoms such as helium or hydrogen migration in irradiated materials.

For the Fe-He system, modeling of helium clusters has been performed by Morishita, Wirth and co-workers (Morishita, Sugano et al. 2003a; Morishita, Sugano et al. 2003b) and Bringa, Wirth et al. (2003)using semi-empirical potentials. In one paper, Morishita, Sugano et al. (2003a) performed molecular dynamics (MD) calculations to evaluate the thermal stability of helium-vacancy clusters in Fe In another paper, Morishita, Sugano et al. (2003b)have looked at dissolution of helium-vacancy clusters as a function of temperature increase using the empirical potentials for the Fe-He system. Wirth and Bringa (Wirth and Bringa 2004) have simulated the motion of one single 2He-3Vac cluster at 1000 K using the same potential system. First principles calculations of helium atoms in interstitial and substitutional sites has been performed by Fu and Willaime 2005 and can be used to provide an input parameter set for kinetic Monte Carlo calculations.

The kinetic Monte Carlo model(Deo, Srivilliputhur et al. 2006; Deo, Okuniewski et al. 2007; Deo, Srinivasan et al. 2007) consists of helium interstitials on the octahedral sublattice and vacancies on the bcc iron lattice. The migration of the free (not clustered) helium is shown in figure 5(a) and that of the vacancies and self interstitial atoms is shown in Figure 5(b). The
lowest-energy migration path of the SIA corresponds to a nearest-neighbor translation- rotation jump in the <111> direction. The rates of migration of the point defect entities are calculated as

i Л

rmigration V migrationexp

where the superscript i refers to the helium, self interstitial atoms and the vacancy point defect entities. The rate of migration of the point defect entity is rmigration, the attempt frequency is v’migration, the migration barrier is Emigration, while kB and T are the Boltzmann constant and the temperature respectively. Two point defect entities are considered to be in a cluster when the distance between them is less than a0, which is the lattice constant of bcc iron. Interstitial atoms are not allowed to dissociate from clusters. Dissociation of the helium and the vacancy from the cluster is described in Figures 2(c) and 2(d) respectively. The rate
of dissociation of a point defect entity (i = helium or vacancy) from a cluster into the bulk lattice is considered to be thermally activated and is calculated as:

( E Л

ri = v exp dissociation (4)

fdissociation udissociationг І і г-,-, I ^ ‘

L kBT J

where Гdissociation is the rate or dissociation, V’dissocation is the attempt frequencies, Edissociation is the energy of dissociation. The dissociation energy Edissociation of a point defect from a cluster

is taken to be the sum of the energy to bind a point defect entity to the cluster and E‘migration.

Small bubbles migrate with a Arrhenius migration rate parameterized by Table 1. Larger bubbles migrate by surface diffusion at the bubble-matrix interface (Cottrell 2003; Evans 2004),

D, = Ds(f^), (5)

where DB and Ds are the bubble and surface diffusivities respectively, Q is the atomic volume and r the radius of the bubble.

Morishita et al. (Morishita et al. 2003) and Fu et al (Fu and Willaime 2005) have calculated the migration energies of helium and vacancies as well as the binding energies of some helium — vacancy clusters. We employ these migration barriers to calculate the rate of migration of the point defect entities. Parameters used to calculate these quantities are described in Table 1. These parameters are used to calculate the rates migration and dissociation events (Equations 1-3) in the system and build the event catalog for the kMC simulation.

The event catalog is generated by calculating the rates of migration or dissociation of the point defect entities using Equations (3-6). The kMC event catalog consists of the migration, clustering and dissociation of the point defect entities, helium, self interstitial atoms and vacancies. The transition probability of each event is proportional to the rate of event occurrence, calculated by the Equations (3-6). We follow the well established kMC simulation algorithm (Bortz, Kalos et al. 1975; Fichthorn and Weinberg 1991) which is a stochastic, atomic — scale method to simulate the time-evolution of defects and nano/ microstructural evolution that focuses on individual defects and not on atomic vibrations.

Reaction pathways available in the system are tabulated in table 1. Rates of migration of events are calculated at each kMC step. Parameters are obtained both from literature using first principles calculations (Fu, Willaime et al. 2004; Fu and Willaime 2005) and also from molecular statics calculations using semi-empirical potentials as employed by Morishita et al (Morishita et al. 2003) and Wirth and Bringa (Wirth and Bringa 2004) (using the Ackland Finnis-Sinclair potential, the Wilson-Johnson potential and the Ziegler-Biersack-Littmark — Beck potential for describing the interactions of Fe-Fe, Fe-He and He-He, respectively). Helium atoms are introduced at random positions at the beginning of the simulation.

Self interstitial atoms (SIA) are produced in the simulation along with vacancies as Frenkel pairs as cascade debris. Self interstitial atoms are mobile and cluster. Self interstitial clusters up to size 5 migrate one dimensionally. Self interstitials of higher size are stationary. Mono — and clustered vacancies are mobile. Vacancies can dissociate from a vacancy cluster as well as a helium-vacancy cluster. Helium atoms migrate on the octahedral sublattice as well as part of helium-vacancy bubbles. A substitutional helium is considered as a 1-1 He1V1 and is mobile. If a bubble has a helium-vacancy ratio greater than 5, it emits a self interstitial atom. Small bubbles migrate according to Equation 3 and Table 1 while large bubbles migrate according to Eq. 3. Self interstitial atoms and vacancies annihilate when they meet either as point defects or in a cluster. The boundary of the simulation cell acts as a sink for the point defect entities.

At each kMC step, the system is monitored to identify a clustering event. When two point defect entities (helium-helium, vacancy-vacancy, helium-vacancy) are in a cluster the simulation creates a mapping between the entities and the cluster such that for each cluster there are at least two entities associated with the cluster. The event catalog is updated with the new rates of event occurrence and the transition probabilities for the next kMC event are calculated using Equations (4-6) using the parameters from Table 1.

Simulations were performed for a damage energy of 100keV. A production rate of randomly distributed cascades is assumed such that the damage is introduced at a rate of 10-6 dpa/ s. The simulation cell size is 400aox400a0x400a0/ where a0 is the lattice parameter of iron. The number of Frenkel pairs introduced in the simulation cell are calculated by the Norgett- Robinson-Torrens relationship,

displacements per cascade = (6)

where Ed is the energy of incident particles while Ed is the threshold energy (40eV for iron (1994)). Initial concentration of helium is a parameter in the simulation and is varied from 1 to 25 appm/ dpa. Damage measured in dpa is another parameter in the system and determines the length of the simulation. The incident energy is also a simulation control parameter that determines the number of defects introduced in a single cascade. The simulation is performed until sufficient damage is accumulated in the simulation cell.

Simulations are performed for values of the He concentrations varying from 1 to 25 appm/ dpa (appm = atomic parts per million). Initial concentrations of both point defect entities are calculated from the NRT formula (Equation 4) using an incident energy of 100keV. Overall damage is introduced at the rate of 10-6 dpa/ s and the simulation is carried out until the required amount of damage has accumulated in the system.

Figure 3 shows the concentration of bubbles as a function of He concentration (appm/ dpa) after damage equal to 0.1 and 1 dpa is introduced. The simulation temperature is 0.3Tm. The bubble density increases with increasing He/dpa ratio. The bubble density can be described by a power law expression,

CB = K (cHe ) , CB = K (cHe) (7)

where is the bubble density, is the helium concentration expressed as appm He/dpa and K and m are constants determined by the kinetic Monte Carlo simulations. We find that the exponent m is approximately 0.5. In Figure 6, the solid lines are fit to the data assuming a square root dependence of the bubble density on helium concentration. Thus the bubble density increases as the square root of the He/dpa ratio. While experimental evidence of this variation is difficult to find, such dependence has been suggested by rate theory calculations (SINGH and TRINKAUS 1992) for the case of cold helium implantation annealing. The square root dependence is also found for equilibrium bubbles containing an ideal gas(Marksworth 1973).

The bubble density increases with damage (expressed as dpa) as damage is accumulated till 1 dpa. At higher dpa ratios more vacancies are produced that may serve as nucleation sites for embryonic bubbles by trapping helium. Thus the bubble density scales directly with increasing damage. That bubble density increases with accumulating dpa is observed in experiments (Sencer, Bond et al 2001; Sencer, Garner et al. 2002; Zinkle, Hashimoto et al. 2003). These experiments suggest much smaller values of bubble density than those calculated in the present simulations. Embryonic bubbles are submicroscopic and difficult to estimate from experimental observations; while the kMC simulations have a large fraction of nanometer size bubbles; making comparison with experiment difficult.

This example demonstrates the use of the kinetic Monte Carlo method to simulate defect diffusion in irradiated materials. The method depends on the ability to calculate rates of migration events of individual defects. The advantage of the method is that all atoms do not necessarily need to be simulated and large time scales are accessible to the simulation provided the underlying physics remains invariant.