Kinetic Monte Carlo simulations of boron activation in implanted Si under laser thermal annealing

The interactions between defects and implanted impurity atoms in semiconductor substrates play a central role in the dopant activation mechanism during annealing processes. In past decades, the study of uniform temperature treatments revealed the relationship between defect kinetics, impurity diffusion, and activation in Si, leading to the discovery of defect-driven transient phenomena.¹⁾ Modified thermal processes are currently needed for ultrashallow junctions in micro- and nanodevices. In particular, submicrosecond laser thermal annealing (LTA), which supplies very localized thermal budgets in a submicrosecond time interval, could satisfy new requirements for junction engineering. Indeed, the localized and transient annealing can lead to activation without diffusion in the solid phase irrespective of whether melting occurs.^2–4) In the molten region, dopant activation is a consequence of impurity trapping,⁵⁾ whereas in the solid phase, it is expected to be the result of ultrafast kinetics in the highly damaged matrix.⁶⁾

The kinetic Monte Carlo (KMC) method should be the key formalism for accurate computational study of this thermally activated kinetics^7,8) in suitably large systems. However, the conventional KMC algorithms developed so far reliably predict the dynamics of dopant-defect system only for the uniform temperature case; consequently, they cannot be applied to the simulation of LTA processes, which induce highly non-uniform thermal fields (the temperature typically drops by a hundred degrees from the hottest regions to the coldest ones) as well as phase changes when melting takes place.

In this paper, we use an original phase-field KMC (PF-KMC) methodology to investigate the kinetics of B implants in Si in the extremely far-from-equilibrium conditions caused by a sequence of implantation and laser irradiation. The particles involved in the system evolution are boron atoms, implantation damage, and boron-interstitial clusters (BICs). A key feature of our method is coupling of the KMC code with a phase-field model,^2,4) which also properly simulates the thermal field evolution in irradiated samples when the liquid–solid phase transition occurs. The PF-KMC simulations correctly predict the experimentally measured local dopant activation levels, demonstrating the central role played by the highly damaged crystalline matrix and the importance of the competitive dopant-defect kinetics in the dopant activation dynamics.

The PF-KMC code applied in this research extends the capabilities of a previous version that was restricted to simulation of only the implantation damage, i.e., interstitial- (I) and vacancy (V)-type defects.⁹⁾ The advantage of the PF-KMC code over conventional KMC schemes is the proper statistical simulation of a system of interacting particles in non-uniform thermal and phase fields T(t, r), Φ(tr), respectively, whereas the conventional KMC works only for a uniform temperature T(t) and homogeneous phase. We note that for the LTA processes considered here (inducing a temperature gradient on the order of ∼1 K/nm), any approximation of T(t, r) using an ad hoc uniform field T(t)^* is not feasible because it dramatically affects the position-dependent rates of the different phenomena considered in the simulation. Hence, it is mandatory to couple the stochastic method with a phase-field model that properly simulates the rapidly varying non-uniform thermal and phase (if melting occurs) fields.⁹⁾

The current version includes dopant atoms that can reside in a substitutional position B_s or form a complex with interstitial defects, i.e., BICs.¹⁰⁾ Without loss of generality, we assumed that I and V point defects as well as the BI complex (one dopant atom bound to one I point defect) are the only mobile species, and that defect clusters X_m (X = I or V) as well as the dopant-defect cluster B_nI_m (a cluster formed by n dopant atoms and m interstitials) can only absorb or emit mobile defects.¹¹⁾ Hence, we introduce the following reactions in our PF-KMC code:

$$\begin{align} \mathrm{I} + \mathrm{B}_{\text{s}} &\leftrightarrow \text{BI}, \end{align} \tag{ 1 }$$

$$\begin{align} \mathrm{I} + \mathrm{B}_{n}\mathrm{I}_{m} &\leftrightarrow \mathrm{B}_{n}\mathrm{I}_{m + 1}, \end{align} \tag{ 2 }$$

$$\begin{align} \mathrm{V} + \mathrm{B}_{n}\mathrm{I}_{m} &\leftrightarrow \mathrm{B}_{n}\mathrm{I}_{m - 1}, \end{align} \tag{ 3 }$$

$$\begin{align} \text{BI} + \mathrm{B}_{n}\mathrm{I}_{m} &\leftrightarrow \mathrm{B}_{n + 1}\mathrm{I}_{m + 1}. \end{align} \tag{ 4 }$$

Concerning the B_n complexes, they can eventually form by Eqs. (2) or (4) and annihilate only through interaction with the mobile BI or I particles.

The free defect migration energies were taken from the work of Pelaz et al.,⁷⁾ whereas the frequency prefactors were fixed to achieve the proper diffusivity values near the silicon melting point ($D_{\text{I}}^{\text{melt}} = 3.6 \times 10^{ - 4}$ and $D_{\text{V}}^{\text{melt}} = 4.0 \times 10^{ - 5}$ cm²/s).¹²⁾ This strategy allows a reasonable physical setting of the free defect parameters in the wide temperature range spanned during the simulated processes. This particular issue makes the LTA simulations more ambitious than those of conventional thermal processes, where the physical properties have to be calibrated in a limited temperature range. The interstitial and vacancy cluster parameters (binding energies and capture radii) are taken from Refs. 11, 13–15. The I_m (V_m) clusters can emit only an I-type (a V-type) point defect, whereas a B_nI_m aggregate can emit a mobile BI or a point defect. The binding energies E_b for the release of a mobile defect from an immobile agglomerate B_nI_m were taken from Ref. 10. In the liquid state, free and clustered defects are assumed to be entirely dissolved, whereas dopant atoms are positioned in a substitutional state B_s during solid regrowth after melting.

To support the theoretical results obtained in the investigation of the boron activation dynamics under LTA treatment, we experimentally studied the B redistribution and its activation in Si under a rather wide range of conditions. Laser irradiation was applied to a boron-implanted silicon substrate (B 30 keV energy, 1 × 10¹⁵ cm⁻² dose, projected range: R_p ∼ 120 nm). The samples were annealed using the EXCICO LTA system (λ = 308 nm) at room temperature, with a pulse duration of 160 ns and frequency of 1 Hz in the multishot configuration. Laser energy densities of 1.5, 2.3, and 2.6 J/cm² were considered for single and multiple pulses (1, 3, 5, and 10 pulses) to achieve different thermal budgets and maximum melt depths R_melt: non-melt R_melt = 0 (1.5 J/cm²), partial melt with R_melt = 75 nm ≤ R_p (2.3 J/cm²), and partial melt with R_melt = 140 nm ≥ R_p (2.6 J/cm²). The melt depths were extracted from the chemical analyses considering the relevant profile modification in the molten region.⁶⁾ After implantation and/or LTA, all samples remained mono-crystalline, as verified by transmission electron microscopy (TEM).

All these processes were also simulated by applying the developed PF-KMC code. The periodic and boundary conditions were set as in Ref. 9. Considering the values of R_melt and R_p, the choice of a simulation box (X_b × Y_b × Z_b) with a depth of Z_b = 600 nm (in the direction of laser irradiation) reliably approximates a bulk sample. The remaining dimensions, X_b = Y_b = 90 nm, were fixed so as to obtain good density resolution in the simulated profiles.

The initial state (post-implant) for the dopant-defect system (i.e., I_m, V_m, and B_nI_m) was derived by the binary collision approximation (BCA) and atomistic KMC simulations of their evolution at room temperature. The simulation results suggest a reduction of 60% in the BCA-estimated total implantation damage stored in small defect aggregates.¹⁶⁾ We note that the initialization of the dopant-defect system for the LTA simulation is a crucial issue. In conventional simulations of dopant-defect systems, it is a common practice to set all the implanted impurity atoms in an interstitial position and, concurrently, set a distribution of self-interstitial atoms equal to the implanted dopant distribution and stored as I₂ clusters.^8,11) It is clear that this choice, which neglects the first stage of the annealing treatment, cannot be implemented for LTA simulations, where the post-implant evolution occurs for only hundreds of nanoseconds.

During the evolution promoted by the laser-induced thermal budget, boron atoms can reach substitutional lattice sites as well as form immobile aggregates with interstitial defects (i.e., BICs). The first component will increase the carrier density, whereas the second will need an additional thermal budget to dissolve and then become active. Figure 1 shows the boron component stored in the B_nI_m clusters for all the cited laser irradiation conditions. Comparing the regions (z > 75 nm) where the samples remain solid for both the 1.5 and 2.3 J/cm² cases, we can observe a pronounced decrease in the local density of BICs with increasing laser fluence. On the other side, a comparable inactive impurity concentration is predicted by KMC simulations when the laser fluence increases from 2.3 to 2.6 J/cm². We demonstrate in the following that this is connected with the similar local activation efficiency of the two processes, denoting the limit of the laser-induced thermal field to release further B atoms from the BICs. Hence, the formed boron-interstitial aggregates play a central role in the solid-phase activation efficiency during the laser-induced kinetics in the highly damaged crystalline matrix.

Zoom In Zoom Out Reset image size

A full atomistic simulation approach for LTA treatment, considering the submicrosecond scale of the process, could be a key pathway for understanding the dopant-defect dynamics through only a detailed comparison with the experimental counterpart. From this point of view, the chemical B profiles were analyzed using secondary ion mass spectroscopy (SIMS), and the electrical activation was evaluated experimentally on both a qualified SSM150 spreading resistance profiling (SRP) system and electrochemical capacitance–voltage (ECV) equipment. In Fig. 2, the SRP (triangles) and ECV (squares) electrically active distributions after LTA are presented and compared to the associated SIMS profiles. Comparisons of the chemical and active (SRP or ECV) profiles demonstrate that significant activation levels are reached under all the annealing conditions, with a relevant exception for the 1.5 J/cm² non-melting case, where activation it is limited to 78% (ECV) and 12% (SRP) with respect to the total implanted boron. The discrepancy between the activation levels extracted using the ECV and SRP techniques can be understood considering the measurement peculiarities. ECV gives direct access to the total charges by integrating with the CV scan the density of states associated with the impurity levels in a wide energy range of the Si gap; the SRP technique involves only the mobile charges in the measurement of the active concentration. The presence of the highly defective matrix in the non-melting process can explain the difference between the ECV and SRP results in terms of both reduced mobility or incomplete activation of the B_s component.

Zoom In Zoom Out Reset image size

Figure 3 shows the absolute value of the percentage difference between the integrals of the solid-phase ECV and SRP curves (i.e., calculated in the region that does not melt) normalized with respect to the ECV integral, again evaluated in the solid region (dark squares and solid line, left scale), whereas the concentration curves are shown as functions of the depth in Fig. 2. The activation levels measured with the two techniques approach each other with increasing laser density energy (i.e., for larger melt depths) when the boron evolution occurs mainly in less damaged regions. The trend of the discrepancy suggests both a reduced mobility and then incomplete release of their own carriers from B_s atoms in the valence band by thermal excitation at room temperature when the matrix is strongly damaged. The red squares and dashed line show the integral of the total solid-phase defects (i.e., I_m, V_m, and B_nI_m) as a function of the laser fluence extracted from the PF-KMC simulations. A strong correlation appears between the defect reduction and the decrease in the discrepancy between the ECV and SRP data. Furthermore, TEM measurements of the same specimens confirmed the presence of dislocation loops, as correctly predicted by the Monte Carlo simulations.

Zoom In Zoom Out Reset image size

The empty squares in Fig. 2 represent the simulated substitutional boron atoms B_s (red). As a reference, the SIMS profiles (dark line) are also shown. The KMC simulation reproduces the full dopant activation in the molten region as well as the measured activation efficiency in the solid phase. The total boron activation efficiency increases from 78% (1.5 J/cm²) to 92% (2.3 J/cm²) and 94% (2.6 J/cm²). In these conditions, the remaining boron atoms are trapped in B_nI_m aggregates. This suggests that the calibration used here, which was conducted for uniform-temperature annealing, correctly predicts the B_s fraction induced by the LTA processes. The activation process induced by this thermal treatment is the result of competitive events, namely, dissolution/capture of the mobile BI from/by a B_nI_m complex and its diffusion in the highly damaged crystalline matrix before it is promoted in a substitutional site by I release. From this point of view, it is fundamental to simulate the full dopant-defect system during LTA in order to predict the correct activation level.

The PF-KMC code also reproduces the experimentally observed rapid activation saturation with increasing shot number. Figure 2 (2.3 J/cm² panel) shows both the experimental ECV data and KMC results. The red and blue symbols represent the substitutional boron profiles after a single irradiation and a 5-pulse process, respectively. To strengthen our analysis, we also measured and simulated, for the three fluence cases, the activation profiles for 3- and 10-pulse processes. The strong overlap between all the experimental and simulated profiles irrespective of the pulse number suggests that the low thermal budget repeatedly supplied to the substrate is not sufficient to effectively release impurity atoms from the more stable B_nI_m complexes formed after the first laser pulse. By comparing the ECV and chemical data under these conditions (2.3 J/cm²), solid-phase activation efficiencies of 74 and 76% was found experimentally for single and 5-pulse irradiations, respectively (the integral is restricted to only the non-molten region). The KMC simulations for the same processes give levels of 72 and 76%, respectively. Hence, a reliable agreement can be noticed also in multipulse annealing, where the activation efficiencies are almost completely saturated after more than 5 laser pulses (not shown).

In conclusion, our PF-KMC simulations of the boron damage evolution during submicrosecond LTA suggest that the defective system stabilization in the initial stage of the process, composed of I and V clusters as well as BICs, plays a key role in the boron activation process. A fraction of the implanted impurity ions also remained trapped in the BICs when the number of shots, i.e., the thermal budget, was increased. However, good activation levels were reached in both the molten and non-molten regions. In addition, our results demonstrated that, although their inner atomic structure is not considered in our simulation approach, the energetics of the I and V defects in silicon as well as the BICs extracted from constant-temperature processes correctly predicts the activation dynamics in a highly damaged matrix during a submicrosecond laser annealing treatment.

Acknowledgments