Kinetic Monte Carlo simulations for transient thermal fields: Computational methodology and application to the submicrosecond laser processes in implanted silicon

PHYSICAL REVIEW E 86, 036705 (2012) Kinetic Monte Carlo simulations for transient thermal fields: Computational methodology and application to the submicrosecond laser processes in implanted silicon G. Fisicaro,1,*L. Pelaz,2P. Lopez,2and A. La Magna1 1CNR-IMM, Zona Industriale VIII Strada 5, I-95121 Catania, Italy 2Department of Electronics, University of Valladolid, ES-47011 Valladolid, Spain (Received 1 June 2012; published 27 September 2012) Pulsed laser irradiation of damaged solids promotes ultrafast nonequilibrium kinetics, on the submicrosecond scale, leading to microscopic modifications of the material state. Reliable theoretical predictions of this evolution can be achieved only by simulating particle interactions in the presence of large and transient gradients of the thermal field. We propose a kinetic Monte Carlo (KMC) method for the simulation of damaged systems in the extremely far-from-equilibrium conditions caused by the laser irradiation. The reference systems are nonideal crystals containing point defect excesses, an order of magnitude larger than the equilibrium density, due to a preirradiation ion implantation process. The thermal and, eventual, melting problem is solved within the phase-field methodology, and the numerical solutions for the spaceand time-dependent thermal field were then dynamically coupled to the KMC code. The formalism, implementation, and related tests of our computational code are discussed in detail. As an application example we analyze the evolution of the defect system caused by P ion implantation in Si under nanosecond pulsed irradiation. The simulation results suggest a significant annihilation of the implantation damage which can be well controlled by the laser fluence. DOI: 10.1103/PhysRevE.86.036705 PACS number(s): 02.70.−c, 61.72.Cc, 61.43.Bn, 82.20.Wt I. INTRODUCTION The kinetic Monte Carlo (KMC) method is the key formalism for the simulation of evolution of large systems when the kinetics is driven by thermally activated phenomena. Its applications cover different topics: e.g., erosion [1], crystal growth [2], surface modification by plasmas [3,4], cell and macromolecule interactions [5], and so on. The KMC method simulates the kinetics by a stochastic sequence of elementary events (e.g., displacement, chemical reaction, bond formation, dissolution, etc.) leading to a state change of the evolving Monte Carlo particles (e.g., atoms, defects, molecules, cells). A common characteristic of the KMC models developed so far is that the temperature is uniform in the simulated space region [i.e., T(r,t)=Tuni(t) for all rin the simulation box]. This assumption allows the stochastic selection of the temporally ordered chain of events using temperature-dependent probability rules which depend only on the event type. We note that the optimal algorithmic formulation and the practical implementation of the KMC model in simulation codes are strongly based on the one-to-one relationship between the event type and its probability (and frequency) value and that this unique correspondence is not in general valid when Tis not uniform. An interesting application of the KMC approach concerns the simulation of the damage evolution in ion-implanted samples [6,7] for technological applications in the field of microelectronics. Several works have demonstrated the prediction power of this method when it is applied to postimplantation annealing processes at high uniform temperature. However, the KMC approach has not been applied so far in the case when the annealing induces, in the damaged region, large and transient gradients of the thermal field, as in the case of pulsed laser irradiation. *giuseppe.fi[email protected].it Simulations of laser processes are crucial in order to understand the effects of irradiation on material modifications in both melting and solid phases. Many issues are also still poorly understood for the widely investigated case of laser thermal annealing of implanted silicon; for instance, the behavior of common dopant atoms (arsenic, boron, or phosphorus) and their evolution in the melting phase or the dopant activation promoted by laser irradiation processes in solid regions. Due to the intrinsic features of these annealing treatments, i.e., localized irradiation with a spaceand time-dependent (on the nanosecond scale) thermal field, experimental analysis can reveal only specimen postirradiation characteristics. In order to understand in depth how the defect system (i.e., the damage) evolves and reaches its final state, a theoretical work based on full process simulation is crucial and KMC simulation could be the ideal method. Conventional KMC algorithms cannot be applied for the simulation of the laser annealing process since it induces a highly nonuniform Tfield and, in contrast to processes at uniform temperature, the spaceand time-dependent temperature implies spaceand time-dependent event frequencies. An evolution and a change of the thermal field requires an update of the frequencies associated with all the possible events. This particularity requires the introduction of different and more efficient sampling methods. In this paper we propose a KMC model which explains the postimplant kinetics of the defect system in the extremely far-from-the equilibrium conditions caused by laser irradiation in the nonmelting, melting, and partial melting regimes. As presented in Sec. II, it considers defect diffusion, annihilation, and clustering. A key feature of our method is the coupling of the KMC code with a phase-field model which properly simulates the thermal-field evolution in an irradiated sample also when the liquid-solid phase transition takes place. In Sec. III we discuss the coupling methodology and in particular the procedure 036705-1 1539-3755/2012/86(3)/036705(11) ©2012 American Physical Society G. FISICARO, L. PELAZ, P. LOPEZ, AND A. LA MAGNA PHYSICAL REVIEW E 86, 036705 (2012) implemented in order to map in the KMC simulation box the fast-varying nonuniform thermal field T(r,t) derived by means of the thermal-field simulations. In Sec. IV we present a continuous approximation of the KMC model which results in a set of reaction-diffusion partial differential equations (PDEs) controlling the average density fields of the MC particles. Comparisons between the particlelike and the continuous simulations are useful for testing and for a proper discussion of the application limits of the continuous approximation. In order to consider a realistic application, the KMC code has been applied to the laser annealing processes of a phosphorus-implanted silicon substrate (P ions with 200 keV energy and a dose of 1014 cm−2). The evolution of point defects and clusters has been analyzed by varying the laser fluence. Simulation results, discussed in Sec. V, predict a reduction of the implantation damage with increase of the applied thermal budget. These results correctly match the increase of the dopant activation efficiency in the silicon solid phase as the laser fluence builds up. The same scenario has been simulated by means of the continuous model. A comparison of the two approaches will supply a test for the PDE model and the correct parameter relations when passing from the KMC to the continuous equation methodology. II. KMC CODE FOR THE DEFECT SYSTEM IN SILICON Dopant-defect system evolution in implanted silicon under constant temperature processes has been widely investigated in recent decades. Consequently, a great number of codes based on KMC and continuous models have been developed, most of them daily used in academic and industrial research [7]. In this section a brief generic description of the KMC method will be given, focusing on the modifications of the numerical formalism which make possible the kinetic simulation in the case of a transient nonuniform thermal field. Indeed, the transition from the conventional annealing to the laser irradiation process imposes selection solutions usually not implemented in the conventional KMC algorithms. The KMC formalism is consistent with the framework of the transition state theory and, therefore, we assume that all the elementary events occur when an activated intermediate state is achieved with the aid of the thermal fluctuations. The event probability (and consequently the event rate or frequency) is a function of Eact/kBT, where Eact is the energy difference between the activated state and the initial one, kBthe Boltzmann constant, and Tthe temperature. Event rates are also dependent on several other parameters (e.g., capture volumes, entropy contribution, etc.) whose values must be defined apriori.We consider an Arrhenius-like Tdependence of the probabilities. Therefore, the frequency (or probability per unit time) related to the occurrence of a particular event ican be expressed as νi=1 τi =ν0 iexp −Eact i kBT,(1) where ν0 iis a prefactor which formally includes the dependence of νion the other parameters. In our (nonlattice) KMC code the damage in Si is described as a system of mobile point defects, interstitials (I), and vacancies (V), and their immobile aggregates (from small Laser irradiaon Laser irradiaon Liquid Liquid Solid Solid FIG. 1. (Color online) Schema of no-lattice kinetic Monte Carlo environment which considers interstitial (dark balls) and vacancy (white balls) defects in mobile (surrounded by arrows) or clustered configuration. The two panels refer to the melting front evolution induced by a laser irradiation process. Particle annhilation in the molten region [red (dark gray)] and defect ripening process in solid phase [yellow (light gray)] are also represented. clusters to extended defects) residing in a simulation box (Xb×Yb×Zb). Both mobile and immobile particles are located in internal positions r(i) of the three-dimensional simulation box. The system status is described as the set of necessary information on the defects (i.e., defect position, type, and size) which is stored in proper arrays. The regular Si atoms residing in the sites of the crystalline silicon lattice (i.e., undamaged regions) do not participate in the evolution; and, as a consequence, their positions are not stored in our code. A two-dimensional schema of the KMC environment is reported in Fig. 1with interstitial (dark balls) and vacancy (white balls) defects in mobile (surrounded by arrows) or clustered configurations. Without loss of generality, we also assume that Iand Vpoint defects are the only mobile species and defect clusters Xncan only absorb or emit point defects. Therefore, in order to reproduce the defect system evolution as a function of time, a numerical procedure simulating the effects of the following elementary events has been implemented: I(r)→I(r)(2) (Ijumps randomly in order to reproduce diffusive behavior, event rate νIdiff); V(r)→V(r)(3) (Vjumps randomly in order to reproduce diffusive behavior, event rate νVdiff); I+V↔0(4) (pair annihilation, event rate νIV ); X+Xn−1↔Xn(5) (Ior Vaggregate growth and dissolution, event rates νXXn−1 and νXndiss); In+V↔In−1(6) (Vannihilation in Iaggregates, event rate νInVνIn−1Vdiss); Vn+I↔Vn−1(7) (Iannihilation in Vaggregates, event rate νVnIνVn−1Idiss) 036705-2 KINETIC MONTE CARLO SIMULATIONS FOR TRANSIENT ... PHYSICAL REVIEW E 86, 036705 (2012) where rand rare the positions occupied by the defects in the simulation box before and after the random jump. If T is uniform in the whole simulation box, there is a one-to-one correspondence between the event type and its associated rate, independently of the position of the particle in the simulation box. As a consequence we can classify the events in discrete classes {CIdiff,CVdiff,...},where each class is characterized by a single common rate {CIdiff ⇒νIdiff,CVdiff ⇒νVdiff,...}. This is the case of KMC algorithms simulating dopantdefect systems in implanted semiconductors during constant temperature annealing. The parametric dependence on T allows the implementation of an efficient sampling method for the selection of the occurring event from the set of possible ones. If the system at a given time could evolve by one of the next possible events i, the total event rate is R= i νi= C∈classes NCvC,(8) where NCis the number of events belonging to the Cclass at that time. The random sampling process is composed of two steps: the first selects a particular class by means of the discrete probability distribution (NCvC/R) among all the possible classes; the second picks randomly one of the NCevents belonging to the selected class. Hence standard selection algorithms can be implemented based on the class storage for each particular elementary event, and the simulation time is increased by an amount t=1/R. Moreover, the update of the system state and the event distribution in the probability classes is simple, since it requires a local change of the particle type and a well-defined modification in a very limited part of the whole information stored in the {NC}array. These algorithms cannot be applied for a kinetic Monte Carlo simulation of the laser annealing process since the temperature is a nonuniform variable field T(r,t) (represented by means of the color scale in the schematic of Fig. 1). In this case, assuming the same functional dependence (1) of the event rate on T(which is valid in the limit of pointlike events), we obtain a continuous field map for the frequencies vi(T(r,t))which vary with time and space in the simulated samples. Therefore, standard selection algorithms, based on class storage, cannot be implemented since events in the same class have different frequencies (i.e., similar particles involved in the transition experience different temperatures) and class ordering does not result in computational advantages. We have to develop a code which considers the possibility of fieldlike event probabilities expressed by suitable local functions of T(r,t). An evolution and a change of the thermal field require an update of all the frequencies associated with all the possible events. This particularity requires the introduction of a different sampling method. We could simply select the event by directly using the probability distribution vi(T(r,t))/R, where the index iruns over all the Nevent events. However, this trivial selection procedure is not efficient for large systems since it requires an o(Nevent) CPU time per KMC step. The use of a tree storage of vi(T(r,t))can overcome these difficulties, and, for consistency’s sake, we will briefly recapitulate this procedure here. A tree is similar to a linked list in being made up of elements and pointers, but each element or node can be accompanied by more than one pointer. The most common type of tree is the binary tree in which each node has two pointers to other descendant nodes and its element wp=w1+w2is the sum of the values present in the descendant nodes. The root node R=iνiis that at the very beginning of the tree while the nodes at the very end store vi(T(r,t))and are called leaves. During a random selection process, a random number between 0 and 1 allows us to choose one of the two descendants of a parent node according to their own probabilities {w1/wp,w2/wp}. Hence starting from the tree root, a loop allows us to reach a leaf at the end of the tree, and the choosing procedure is equivalent to the direct use of the vi(T(r,t))/R weights. The number of steps needed to complete this selection is proportional to the logarithm (base 2 in the case of a binary tree) of Nevent [o(logNevent)CPU time per step]. Therefore, it is more efficient with respect to the selection based on the simple array storage. The tree root, being the sum of all the event probabilities, i.e., the total event rate Rof Eq. (8), is used to update the time interval between two successive KMC steps. The following procedures are implemented to simulate diffusion and dissolution events. If the selected event is a random jump (i.e., a diffusive event) of the mobile point defects (Ior V), the code updates the particle position from rto r=r+dr,simulating a jump in a random direction. In this case the module of dr[dr=(dr,θ,ϕ) in polar coordinates] is constant and equal to the capture radius of the I-Vannihilation, while θand ϕare randomly picked. If the selected event is the dissolution of an immobile aggregate Xn, the code updates the cluster size from nto n−1 and, concurrently, creates a new mobile point defect at r=r+dr.Here ris the original position of the aggregate and dris again a vector with a random direction and fixed length (equal to the sum of the capture radius of the I-Vannihilation and the capture radius characteristic of the cluster Xn). If the moved particle has neighboring ones to interact with in its new position, a new interaction is added in the possible list. Kinetic Monte Carlo codes based on crystalline lattices can quickly check the present of neighboring particles to interact with, looking at next-neighbor or next-nearest-neighbor sites in the lattice structure [8]. Nonlattice KMC codes are more efficient as they deal with only the defect system, neglecting the local lattice details, but they lose the local correlation between the interacting particles. In principle, a complete checking among the whole set of particles has to be carried out, and for distances greater than their capture radius. In order to overcome this problem, which worsens with increase of the total particle number, a linked list has been implemented, as presented in Ref. [9]. It creates a local correlation between the interacting particles, replacing the lattice structure in lattice KMC codes. We end this section by describing the implementation details both for the interaction mechanisms and for the derivation of the KMC model parameters. A space cutoff is considered for the interactions (4)–(7) (i.e., νi=0 if the particle distance is larger than rn). Therefore, I(V) annihilation, as described by Eqs. (4),(6), and (7), has been assumed when a mobile defect falls, after a random jump, inside the capture radius rn(I−Vn)[rn(V−In)] of another point or clustered V-I-type defect. On the other hand, cluster growth [Eq. (5)] occurs when a mobile point defect falls inside the capture 036705-3 G. FISICARO, L. PELAZ, P. LOPEZ, AND A. LA MAGNA PHYSICAL REVIEW E 86, 036705 (2012) radius rn−1(Vn−1)orrn−1(In−1) of an immobile aggregate of the same kind. Point defect and cluster parameters for both I-type and V-type defects were taken from fundamental theoretical works and experimental literature. In particular, for diffusive events, the activation energy barrier Eact iof Eq. (1) represents the migration energy while the ν0 ifrequency prefactor is related to the D0exponential prefactor of the experimental Arrheniuslike form by means of the relation f=6D0 λ2,(9) where λis the jump distance. Concerning the cluster dissolution events, the activation energy barrier Eact inow represents the dissolution energies, and a probability prefactor equal to the diffusive one (distinguishing between Iand Vdefects) multiplied by the cluster size was chosen. III. COUPLING BETWEEN KMC METHOD AND THE PHASE FIELD The innovation of our code lies in the extension of the KMC method to deal with the evolution problem for the damage in silicon in a spaceand time-dependent thermal field. However, the T(r,t) field must be calculated with an external code in order to simulate the transient conditions of a real process. Pulsed laser irradiation can induce a strongly varying T(r,t) (on the submicrosecond and submicrometer scale); moreover, if the energy density exceeds a proper threshold value, a localized melting can occur. Therefore, the evolution of the molten region (phase field) should also be simulated and reliably mapped in the KMC code. Mapping of the thermal and phase fields induced by the laser annealing processes has been obtained by coupling the KMC code with simulation results carried out within the phase-field methodology. According to this approach a continuous field ϕ(−1⩽ ϕ⩽1, with ϕ=−1 and ϕ=1 in the pure liquid and solid phases, respectively) is introduced to describe smooth phase transitions. The phase-field and heat equations are solved concurrently using the following equations [10]: τ∂ϕ ∂t =W2∇2ϕ−ϕ(ϕ2−1) −λcp Lfus (T−TM)(ϕ2−1)2, (10) ρcp ∂T ∂t −∇ 2(KT)=ρLfus 2 15 8(ϕ2−1)∂ϕ ∂t +S(r,t), (11) where TMand Lfus are the melting temperature and the latent heat of crystalline silicon, respectively. S(r,t) is the heat source due to the absorption of laser light which can eventually be calculated by self-consistently solving the Maxwell equations [11,12]. All the material properties (heat capacity cp, heat conductivity K, density ρ) are phase dependent. The phasefield parameters (τ,W, and λ) are chosen in order to correctly achieve the thin-interface limit [10]. The phase-field model calibration is equivalent to that reported in Ref. [13]. The defect evolution should not substantially influence the T(r,t) and ϕ(r,t) solutions, and therefore the parameters used in Eqs. (10) and (11) do not depend on the damage state. As a consequence, a sequential coupling between the phase-field and KMC models has been implemented as follows. Once a particular laser process has been simulated with the phase-field model, the temperature field as a function of time and melt depth evolution in the processed region is extracted and stored at constant intervals with a given time resolution δt. We note that the minimal value of δtmin is given by the inner discretization of the time variable in the PDE solver of Eqs. (10) and (11) and that the choice δt =δtmin is obligatory in the eventual case of interdependence between the optical and thermal parameters of the phase-field model and the damage state. However, due to the assumed lack of interdependence, a different value δt > δtmin can be used if the defect stage is slowly varying on the time scale δtmin. The choice of δt is related to the time variation of the event frequencies by the local temperature rate δT , Eact min kBT2δT < δR R<Eact max kBT2δT, (12) whereEact min (Eact max) is the smallest (largest) activation energy involved in the event frequencies of the process. Equation (12) has been derived from Eqs. (1) and (8). For instance, for our defect parameter setting we have Eact max =E8 diss =4.22 eV,i.e., the dissociation energy of I8(see Table II). Therefore a temperature variation of δT ∼10 K from the room temperature value leads to a total Rvariation of less than ∼5.4%, while at higher temperature, e.g., T=1600 K, the total Rvariation is less than ∼2.0%. The space dependence of the temperature fields at different times is stored in a grid representing equally spaced points of the simulation box; therefore a linear interpolation allows us to calculate the temperature for the particles occupying a given position, and the related local frequency can be properly updated. Note that during the δt interval, the frequencies are considered only space dependent. Concurrently, by means of the storage of the evolving melt depth curve, the timedependent molten region can be tracked, and all the particles that during their evolution fall inside this zone are fully annihilated. The melting process evolution and the correlated effects in the no-lattice kinetic Monte Carlo environment are represented in the schema of Fig. 1. The two panels refer to the melting front evolution induced by a laser irradiation process. In particular, particle annhilation in the molten region (red) and the defect ripening process in the solid phase (yellow) are represented. Considering the laser pulse shape and analyzing the simulation output, for all the laser processes considered in the following, a value of δt =0.5 ns ensures temperature variations lower than 10 K with a consequent acceptable variation of the total R. IV. FULLY CONTINUOUS MODELING Several advantages arise from the use of the kinetic Monte Carlo approach to simulate dopant-defect system evolution on semiconductor substrates. One of them is the absence of limitations in the number of interactions between dopant atoms and defects. The detailed dopant clustering can be considered 036705-4 KINETIC MONTE CARLO SIMULATIONS FOR TRANSIENT ... PHYSICAL REVIEW E 86, 036705 (2012) as well as dopant-defect cluster complexes without restriction on the number of constituents. The kinetics is simulated by a sequence of “virtual” events associated one to one with “real” events; and, in general, we do not need to reduce the complexity of the problem under study since the computation speed does not slow down with increasing number of interactions considered. At the same time, consistently, the greater part of the parameters can be directly extracted from atomistic (i.e., ab initio) calculations. However, CPU time grows with the total number of MC particles and this fact imposes limitations on the size of the system which can be simulated. Hence, along with KMC methodologies, great interest has developed in the past decades in continuous models, which deal with the same issues previously reported [14–16]. The continuous model, analogous to the KMC method introduced in the previous sections, can be formulated by coupling the phase-field model with a set of rate equations for the defect evolution in the following form [17]: ∂CX ∂t =−∇JX+RX,(13) where RXaccounts for the point defect I-Vand point-defect– cluster interactions [Eqs. (2)–(7)]. In particular, RX=−RI,V +RC,(14) RI,V =kf I,V (CICV−C∗ IC∗ V),(15) kf I,V =4π(DI+DV)rcapt(I,V),(16) where CX,C∗ X,and JXare respectively the density, equilibrium density, and flux of a defect X(Ior V). Possible reactions between Xand other species (point defects, clusters) are taken into account by RX(see Ref. [18] for details). The flux term JX is the sum of the diffusion term (Fick’s law) and a thermal-field drift term [19]: JX=−DX∇CX+CX ∇T T.(17) The diffusivities are written in a general Arrhenius-type expression [20]. Concerning clustered defects, a reduced model, able to capture the main features of the kinetic damage in Si during an irradiation process [16], has been adopted, where clusters with sizes up to 9 are taken into account explicitly with specific dissolution energies: ∂CXn ∂t =kf nCXCXn−1−kb nCXn+kb n+1CXn+1−kf n+1CXCXn. (18) Two different set of reaction rates kf nand kb nhave been tested. The first set [18] represents the conventional choice implemented in most of the continuous models which describe dopant-defect systems under constant temperature annealing: kf n=α4πrnDX,(19) kb n=βDXC∗ Xrn−1exp Q−En diss kbT,(20) where rnand En diss =Em+En binding are the capture radius and the dissolution energy of the cluster Xn,Qthe activation energy of the product DXC∗ X, and αand βcalibration parameters. Emis the migration energy of the mobile defect produced by the dissolution event and En binding is the binding energy of the n-size cluster. The second set, taken from the work of Rafferty et al. (see Ref. [21]), does not include additional constants or the equilibrium point defect concentration C∗ Xin the backward reaction rate kb n: kf n=4πrnDX,(21) kb n=nDX(a)−2exp −En binding kbT,(22) where ais the average interatomic spacing fixed equal to the rcapt value described before, i.e., the capture radius of I-Vannihilation. In both approaches, the bigger clusters XCl are considered to be similar to extended defects with a unique reaction rate, and cluster diffusion, as in Ref. [19], is neglected. It is worthy of emphasis that in the backward reaction rates kb nthe energies in the Arrhenius-type expressions are fixed equal in both sets, while differences appear only in the exponential prefactors. The prefactor for interstitial defects in the case of Giles et al.’s set is ∼2.7×104greater than that of Rafferty et al. For vacancy-type defects the ratio is ∼2.1×102. V. SIMULATION RESULTS AND DISCUSSION A. Calibration and initialization of KMC and PDE models As a test case for study, the damage evolution in implanted silicon for three different laser irradiations has been simulated, reproducing real experimental conditions [16]. In this experiment silicon substrates were implanted with phosphorus ions with an energy of 200 keV and a dose of 1014 cm−2. Subsequently, the specimens were irradiated by means of an excimer laser (308 nm wavelength). Different irradiation processes have been reproduced, varying the laser fluence: with values of 2.6 and 3.0 J/cm2, partial melting processes (melt depth smaller than the P projected range) were obtained. Partial melting processes represent the ideal test case for our KMC code since the implantation damage, residing in the solid region below the maximum melt depth, undergoes the most critical conditions (proximity of a liquid region, large transient gradient of T). On the other hand, at 3.6 J/cm2 a total melting process (melt depth greater than the P projected range) occurs and most of the implantation damage is than fully annihilated. As described in Ref. [16], after implantation and/or laser irradiation (LI), all the samples remain monocrystalline as verified by transmission electron microscopy. From a simulation point of view, this allows implementation in the code of both the physics and the complete set of process parameters for the case of a crystalline silicon substrate. Since the scope of this section is to supply a comparison and a cross validation of the two simulation approaches, without lack of generality, we have tested our KMC code in a calibration setting which is similar to that already used in Ref. [16] for the PDE modeling of the same problem. Hence the same energetic and capture radii have been implemented 036705-5 G. FISICARO, L. PELAZ, P. LOPEZ, AND A. LA MAGNA PHYSICAL REVIEW E 86, 036705 (2012) both in the KMC code and in the continuous model, and the difference arises only in the implementation details of the physical phenomena peculiar to each approach. The reference experimental results refer to nonpatterned (i.e., onedimensional) samples; as a consequence, periodic conditions for the diffusion events and particle interactions were set in the directions orthogonal to the depth. Annihilation conditions were fixed at the irradiated surface to introduce the physical effect of mobile point defect dissolution at the sample surface. Reflective conditions were set on the opposite side (to the surface) in order to reproduce with a good approximation the real bulk behavior. This assumption is reliable when the simulation box length Zbis much larger than the size of the damaged region. The latter region is fixed by the initial conditions as well as the annealing process (temperature and duration) conditions. The projected range of the phosphorus ions (200 keV) is 300 nm and the maximum melt depth corresponding to the highest laser fluence employed (3.6 J/cm2) is 350 nm, so the choice of Zb=1μm ensures the validity of the real bulk approximation. On the other hand, the choice of the remaining dimensions Xband Ybof the simulation box is determined by the total number of particles to simulate in a KMC run. This number, which plays a central role in CPU time, can be determined by the initial defect density, strictly related to the implantation process characteristics. In our case, a simulation box of Xb=50 nm, Yb=50 nm, and Zb=1.0 μm implies that 13 ×106particles must be simulated. Free defect migration energies were taken from the work of Pelaz et al. [6], whilst the D0values, which also determine the frequency prefactors for the random jump events by means of Eq. (9), were fixed in order to achieve the proper diffusivity values near the silicon melting point (Dmelt I=3.6×10−4 and Dmelt V=4.0×10−5cm2/s[22]). This strategy allows a reasonable physical setting of free defect parameters in the whole temperature range reached during laser annealing processes. This last issue makes simulation of systems treated by laser irradiation, from a calibration point of view, more ambitious with respect to conventional temperature annealing where physical properties of the simulated system have to be calibrated in a limited temperature range. Migration energies, D0values, and point defect equilibrium concentrations are reported in Table I. The products DXC∗ Xfor both interstitial and vacancy point defects have been taken from [23], where activation energies of 4.95 eV (I) and 3.80 eV (V)have been extracted. Interstitial cluster parameters (dissociation and binding energies) are taken from Ref. [18], while vacancy cluster parameters were fitted on experimental supersaturation data [24]. The cluster energetics are reported in Table II TABLE I. Parameters used in the kinetic Monte Carlo code and continuous model for Iand Vpoint defects [6,22,23]. Expression C∗ I43.0×1025 exp( −3.95 eV/kT )(cm −3) C∗ V33.9×1023 exp( −3.3 eV/kT )(cm −3) DI0.34 exp( −1.00eV/kT )(cm 2/s) DV1.23 ×10−3exp( −0.5 eV/kT )(cm 2/s) rcapt 2.725 ×10−8(cm) TABLE II. Iand Vcluster parameters used in the kinetic Monte Carlo code and continuous model [16,18]. InEdiss (eV) VnEdiss (eV) rn(˚ A) X23.61 1.7 1.92 X33.74 2.00 3.84 X43.92 2.20 5.76 X53.79 1.80 7.68 X63.58 3.46 9.60 X73.59 2.40 11.52 X84.22 2.50 13.44 X93.54 2.80 15.36 XCl 4.13 3.46 15.60 together with the corresponding capture radii implemented both in the KMC code and in the continuous one. In the liquid state, free and clustered defects are assumed to be entirely dissolved in the KMC approach, while in the continuous model free defects are forced to be at their equilibrium concentrations at the melting temperature and clusters are assumed to be fully annihilated. The implementation in the KMC approach of the equilibrium effect (i.e., the defect density balance controlled by the local thermal field) is not a trivial question. Nevertheless its effect can be neglected in our case since the local supersaturation maintained by clusters is orders of magnitude greater than the equilibrium value. The defect system was initialized by coupling binary collision approximation (BCA) and atomistic KMC simulations of the defect evolution at room temperature after the implant. The dark dashed (I-type defects) and dotted (V-type defects) lines of Fig. 2refer to the as-implanted defects obtained by stopping and range of ions in matter (SRIM) simulations [25]. In the same figure, the inset reports their difference, which correctly shows higher vacancy concentration near the implanted surface and a prevalence of interstitial defects in the implant tail. Green (light gray) lines refer to the residual defect density after room temperature evolution obtained by means of the atomistic 0.0 0.1 0.2 0.3 0.4 0.5 10 18 10 19 10 20 10 21 10 22 Density [cm-3] Depth [ µ m] 0.0 0.2 0.4 -3x1019 -2x1019 -1x1019 0 1x1019 2x1019 3x1019 Density [cm-3] Depth [ µ m] V-I SRIM data FIG. 2. (Color online) Total interstitial (dashed) and vacancy (dotted) defect density: dark lines refer to the as-implanted defects obtained by SRIM simulation; green (light gray) and red (dark gray) lines refer, respectively, to the residual defect density after room temperature evolution obtained from the atomistic KMC simulation and from the continuous PDE model. The inset shows the difference between total I-andV-type defects as obtained after SRIM simulation. 036705-6 KINETIC MONTE CARLO SIMULATIONS FOR TRANSIENT ... PHYSICAL REVIEW E 86, 036705 (2012) KMC simulation (i.e., the constant Tcase simulating the effective postimplantation residual damage). The simulation results suggest a reduction of 40% of the BCA-estimated total implantation damage. The postimplantation damage is stored in small defect aggregates [26] corresponding to space-dependent cluster size distributions of Iand V-type aggregates. Note that the Iand Vdistributions are always superimposed due to fast rates of the point-defect–pointdefect or point-defect–cluster interactions with respect to the particle diffusion. Similar postimplantation evolution at room temperature has been carried out by means of the continuous model described in detail in Sec. IV. The red (dark gray) dashed (I-type defects) and dotted (V-type defects) lines of Fig. 2represent the residual defect density after room temperature evolution obtained within this approach. The continuous equations of Sec. IV as well as the model for the thermal problem of Sec. III have been solved within an external PDE solver [27]. Excellent agreement can be noticed with respect to the atomistic KMC data [green (light gray) lines]. The comparison provides also a test for the continuous model applied to the case where I-Vannihilation or cluster aggregation are the main kinetic events, since cluster dissolution cannot take place at room temperature. In the following we turn to the investigation of defect system evolution during a laser annealing treatment where all interactions play crucial roles. In these conditions defect complexes can dissolve because the temperature, in the irradiated specimens, reaches the melting point in silicon solid phase near the liquid-solid interface. B. Irradiation simulations In order to simulate laser irradiation processes with the KMC code, the melt depth and thermal fields at different sampling positions (every 0.1 μm starting from the irradiated surface) have been extracted from the phase-field simulations. Figures 3and 4report typical temperature and phase evolutions during a melting laser process. In particular, the time evolution of the temperature at three different positions for the 2.6 J/cm2 laser process is shown in Fig. 3. Note that, after 130 ns from the beginning of irradiation (see the inset in Fig. 3), the maximum 0.0 0.5 1.0 1.5 2.0 2.5 400 600 800 1000 1200 1400 1600 Temperature [K] Time [µs] 0.0 0.2 0.4 0.6 100 300 500 700 900 T variation [K] Time [µs] FIG. 3. Time evolution of the temperature at three different positions for the 2.6 J/cm2laser process: solid, dashed, and dotted lines refer, respectively, to the thermal field at 0.0, 0.5, and 1.0 μm below the irradiated surface (0.0 μm). The inset shows the evolution of the temperature variation between the points 0.0 and 1.0 μm. 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 Melt depth [ µ m] Time [ µ s] 2.6 J/cm 2 3.0 J/cm 2 3.6 J/cm 2 FIG. 4. Time evolution of the melt depth during the melting laser processes with 2.6 (dotted), 3.0 (dashed), and 3.6 (solid) J/cm2laser fluences. variation of ∼960 K is simulated between the temperature at the surface and that reached at a 1.0 μm depth. Note that these large thermal gradient variations with time have to be taken into account if extended defect complexes (>10 nm wide) characterize the damage evolution. In the latter case, the pointlike approximation, assumed in the present work, fails and an on-lattice formalism has to be implemented [27]. Comparing Figs. 3and 4, it can be seen that the largest rates of temperature change mark the beginning of the melting process. As a consequence these rate values have been considered to set the update time interval δt for the whole irradiation process. For times larger than 500 ns, the specimen is quenched with an almost constant temperature field as evidenced by the overlap of the three curves of Fig. 3.In Fig. 4we show the time evolution of the melt depth during three melting laser processes (2.6, 3.0, and 3.6 J/cm2). These curves indicate the melting process duration and the extent of the molten region where defects (mobile and clustered) simulated in the KMC module are fully annihilated. Considering this analysis of the thermal and phase fields, during the KMC simulations an update of melt depth and temperature is made at each point. A linear interpolation allows for a sufficiently accurate determination of the temperature as a function of the depth in the space regions between two consecutive sampling points. We set an update time of 0.5 ns which maintains the temperature increase at all the sampling points smaller than 10 K for the whole simulated time interval. Figure 5shows the total interstitial and vacancy defect density profiles before (dark) and after [colored (grayscale)] laser irradiations obtained by means of atomistic KMC simulations. Red (dark gray) circles (2.6 J/cm2) and green (light gray) triangles (3.0 J/cm2) refer to two partial melting processes where the corresponding maximum melt depths (108 and 194 nm) are lower than the P ions’ projected range (300 nm). The melting front extent has been extracted from Fig. 4. The code reproduces full damage annihilation in the molten regions due to the liquid-solid epitaxial regrowth, while defect reduction in the solid phase (colored squares) is less pronounced. This result is strictly related to the implemented energetics for migration and dissolution events. The high temperatures reached in the solid region of the sample near the liquid-solid interface and the low point defect supersaturation 036705-7 G. FISICARO, L. PELAZ, P. LOPEZ, AND A. LA MAGNA PHYSICAL REVIEW E 86, 036705 (2012) 0.0 0.1 0.2 0.3 0.4 0.5 10 17 10 19 10 21 Depth [ µ m] Defect density [cm-3] 10 17 10 19 10 21 FIG. 5. (Color online) Postanneal total interstitial (filled symbols) and vacancy (empty symbols) defect densities for 2.6 [red (dark gray)] and 3.0 [green (light gray)] J/cm2laser fluences. Dark squares refer to the simulated total interstitial and vacancy defects after KMC room temperature evolution. Dashed (I) and dotted (V) lines refer to the continuous simulation carried out by means of the Giles et al. (upper panel) and Rafferty et al. (lower panel) parameters for cluster dissolution rates. make cluster dissolution events the most relevant interactions. Subsequently, the point defects detached from the clusters can migrate and fall inside the capture radius of composite defects with the consequent preservation or decrease of the total defect density if the capture results in cluster growth or in I-Vannihilation. The atomistic KMC simulations of Fig. 5 suggest that interactions which preserve the total defect system density prevail over the others. In order to investigate these effects, a deeper analysis has to be carried out which looks inside the defect system. Therefore we chose to follow the evolution of small size aggregates since their lower dissolution energies (see Table II) make them the most relevant particles for the characterization of the system evolution. Analyzing the interstitialdefect case of Figs. 6 and 8, atomistic kinetic Monte Carlo simulations predict, 10 15 10 17 10 19 10 21 I Defect density [cm-3] 0.0 0.1 0.2 0.3 0.4 0.5 10 15 10 17 10 19 10 21 Depth [ µ m] FIG. 6. (Color online) Interstitial concentrations after room temperature (empty squares) and one-pulse (filled squares) laser irradiation of 2.6 J/cm2laser fluence obtained by means of the KMC code. Dashed lines refer to PDE simulation for the same process carried out with the Giles et al. (upper panel) and Rafferty et al. (lower panel) parameters for the cluster dissolution rates. Dark, red (dark gray), and green (light gray) lines refer, respectively, to distributions of clusters formed by two, three, and four defects. in the whole simulated solid phase zone, a weak reduction of the small cluster densities for the two laser irradiations considered in the partial melting regime. Note that different colored (grayscale) squares refer to cluster densities formed by two, three, and four particles, respectively. Concurrently, vacancy aggregates undergo a pronounced Oswald ripening process with a reduction of cluster density with respect to the postimplant concentrations for V2and V3and a relative shift to bigger size complexes. These results suggest a faster dissolution of V-type defects with respect to I-type ones, reflecting the dissolution energy differences for interstitials and vacancies (Table II). The more rapid healing of the damage stored in Vaggregates has a double effect: it decreases the population of small V-type clusters and increases the population of small I-type ones due to the V-Ininteractions. This picture of the damage system evolution can be related to the dopant activation in the silicon solid phase during a laser annealing treatment. As reported in detail in previous work [16], the solid phase dopant activation efficiency grows with increase of the laser fluence (i.e., the maximum melt depth). This activation increase has been related to the damage reduction caused by thermal annealing, which can trap mobile dopant particles during its evolution. Atomistic KMC results, elucidating the damage reduction mechanism, suggest that small cluster dissolution plays a central role in dopant activation, decreasing the trap sites controlled by the Oswald ripening process. Another aim of this research is to analyze how a suitable calibration of the continuous models can reproduce reliable atomistic KMC results. The lines of Fig. 5show post-laserirradiation total interstitial and vacancy defect density profiles. The continuous model results reported in the upper panel were obtained with the Giles et al. parameters [18]forthe cluster reaction rates, while the simulation outcomes in the lower panel are obtained with the calibration suggested by Rafferty et al. in Ref. [21]. A substantial difference can be seen between the continuous model (dash-dotted lines) and the 10 15 10 17 10 19 10 21 V Defect density [cm-3] 0.0 0.1 0.2 0.3 0.4 0.5 10 15 10 17 10 19 10 21 Depth [ µ m] FIG. 7. (Color online) Vacancy concentrations after room temperature (empty squares) and one-pulse (filled squares) laser irradiation of 2.6 J/cm2laser fluence obtained by means of the KMC code. Dotted lines refer to PDE simulation for the same process carried out with the Giles et al. (upper panel) and Rafferty et al. (lower panel) parameters for the cluster dissolution rates. Dark, red (dark gray), and green (light gray) lines refer, respectively, to distributions of clusters formed by two, three, and four defects. 036705-8 KINETIC MONTE CARLO SIMULATIONS FOR TRANSIENT ... PHYSICAL REVIEW E 86, 036705 (2012) 10 15 10 17 10 19 10 21 I Defect density [cm-3] 0.0 0.1 0.2 0.3 0.4 0.5 10 15 10 17 10 19 10 21 Depth [ µ m] FIG. 8. (Color online) Interstitial concentrations after room temperature (empty squares) and one-pulse (filled squares) laser irradiation of 3.0 J/cm2laser fluence obtained by means of the KMC code. Dashed lines refer to PDE simulation for the same process carried out with the Giles et al. (upper panel) and Rafferty et al. (lower panel) parameters for the cluster dissolution rates. Dark, red (dark gray), and green (light gray) lines refer, respectively, to distributions of clusters formed by two, three, and four defects. KMC results (filled and empty symbols) in the upper panel. Indeed, in spite of the excellent agreement between the KMC and PDE models calibrated with the Giles et al. parameters for the case of damage reduction after the implant process at room temperature in Fig. 2. The KMC simulations predict a lower damage reduction with respect to the PDE model for the total Ior Vdensity after one-pulse laser irradiation for all the fluences considered. On the other hand, the Rafferty et al. parameters lead to an excellent agreement between the PDE and KMC total damage annihilation after laser irradiation. In Figs. 6–9comparisons between the interstitial and vacancy small cluster densities after one-pulse laser irradiation are reported. As in Fig. 5, the continuous model results reported in the upper (lower) panels were obtained with the Giles et al. (Rafferty et al.) setting procedure for the cluster reaction 10 15 10 17 10 19 10 21 V Defect density [cm-3] 0.0 0.1 0.2 0.3 0.4 0.5 10 15 10 17 10 19 10 21 Depth [ µ m] FIG. 9. (Color online) Vacancy concentrations after room temperature (empty squares) and one-pulse (filled squares) laser irradiation of 3.0 J/cm2laser fluence obtained by means of the KMC code. Dotted lines refer to PDE simulation for the same process carried out with the Giles et al. (upper panel) and Rafferty et al. (lower panel) parameters for the cluster dissolution rates. Dark, red (dark gray), and green (light gray) lines refer, respectively, to distributions of clusters formed by two, three, and four defects. rates. The analysis of distributions of clusters formed by two, three, and four particles leads to similar conclusions for the calibration comparison. Implementing the Rafferty et al. reaction rates, and in particular the backward reaction rates, an excellent agreement is found between PDE and KMC results for interstitial distributions (lower panels). Going to the V-type distributions, Rafferty et al.’s setting improves the general trend with less pronounced damage annihilation with respect to the Giles et al. setting. It is worth underlining the differences between the Giles et al. and Rafferty et al. calibration strategies. In both approaches, the energies in the Arrhenius-type expressions for the backward reaction rates kb nare fixed equal, while differences appear only in the exponential prefactors. The prefactor for interstitial defects in the case of Giles et al.’s parameter set is ∼2.7×104larger than the Rafferty et al. value. For vacancy-type defects the ratio is ∼2.1×102.This difference is at the origin of the discrepancies between the approaches. From this point of view, the present analysis marks an important step in the development of continuous models for the simulation of nonequilibrium evolution of systems which undergo laser annealing processes on the nanosecond scale: a formulation which fixes reaction rates, forcing a fast equilibration of the defect density (as in the Giles et al.’s case), correctly reproduces the damage for the case of conventional thermal annealing processes while it fails for the simulation of laser irradiation processes, whose behavior can be accurately investigated by means of the atomistic KMC approach. The continuous results also suggests that a fine tuning of the Rafferty et al. rates can more easily match the KMC predictions. VI. CONCLUSION In the present work a kinetic Monte Carlo method has been proposed which simulates the evolution of the damage caused by an ion implantation process in silicon during a laser annealing treatment. Thermal fields were previously obtained and extracted from a simulator based on the phase-field methodology. The atomistic KMC simulations suggest that interactions which preserve the total defect system prevail. The results suggest a faster dissolution of V-type with respect to I-type defects, reflecting the dissolution energy differences between interstitial and vacancy complexes. In particular, small I-type aggregates undergo a negligible solid phase reduction for both partial melting laser irradiations tested, while vacancy complexes show a pronounced Oswald ripening process with a reduction of cluster density with respect to the postimplant concentrations for V2and V3and a relative shift to larger complexes. The picture of the damage system evolution has been related to the dopant activation in the silicon solid phase during laser annealing [16]. Atomistic KMC results, elucidating the damage reduction mechanism, suggest that small cluster dissolution plays a central role in dopant activation, decreasing the number of trap sites controlled by the Oswald ripening process. Concurrently, a comparison between the KMC and continuous models with the two different sets of reaction rates (taken, respectively, from the works of Giles et al. [18] and Rafferty et al. [21]) has been reported which suggests that 036705-9