A Comprehensive Solution for Simulating Ultra-Shallow Junctions: From High Dose/Low Energy Implant to Diffusion Annealing



The persistent semiconductor technology trend of shrinking down device size requires development of veryaggressive technological setups consisting in high dose/low energy implants, followed by rapid thermal anneals (RTA). Since it is now well established that phenomena involved in ion implantation will play a key role intransient enhanced diffusion (TED), advanced process simulation tools need to model accurately these two steps for deep submicron CMOS technologies. This paper presents the comprehensive solution fully integrated into the commercial SILVACO TCAD suite needed by advanced process technology designers. The implantation step simulated with a Monte-Carlo code based on the Binary Collision Approximation (BCA), provides initial impurity and defects profiles. The diffusion simulation includes the latest physical model developments: dopantdefect pair, interstitials or vacancy clusters and mixed dopant-defect clusters formation and evolution. Although the complete model covers a wide range of physical interactions, it is straightforwardly usable thanks to a comprehensive set of parameters obtained using an original calibration methodology.


Implant and Diffusion Models

Ion implantation and radiation damage are modelled by means of a simulation technique based on the BCA. The principal assumption of this approximation is that the interaction of energetic particles may be separated into a series of two-body encounters. The real benefit of this approach is the moderate speed of calculation combined with the possibility of including single crystal structures in the calculation. The slowing down of energetic particles is a result of nuclear and electronic stopping. An universal Ziegler-Biersack-Littmark (ZBL) screening function is used for the interatomic potentials while the electron stopping is modelled through local and non-local contributions to the inelastic energy losses. The non-local is that of Wang, et al. [1], and the local stopping is calculated using the technique of Azziz, et al. [2], with a correction for high energies when energy transfer diminishes as described in [3].

In order to simulate the radiation induced damage the model accounts for all collision events within the collision cascade. As demonstrated elsewhere [4], up to 80% of the created I-V pairs can annihilate due to recombination, so the local arrangement of the displaced atoms and vacancies, from where these atoms originated, becomes important when spontaneous recombination has to be taken into consideration. Furthermore, it is not only the number of defects but also their initial spatial distribution, which influences the residual damage. This pertains especially to rearrangement of defects (formation of clusters, loops, voids, etc.) before subsequent migration and thermal annealing is precluded. As a first approach to the initial distribution of defects, at the completion of the cascade, the model of Snyder and Neufeld [5] is used. The value of the vacancy capture radius was carefully chosen so that the amount of residual damage is comparable to that calculated with other techniques, e.g. Kinetic Monte-Carlo.

The full advanced diffusion model will use as initial input, the impurity and defects profiles calculated previously with the BCA code. This model is made of three separated parts, each one corresponding to one particular phenomenon, dynamically interacting with each other [6].

The Classical Dopant Diffusion model [7] takes into account all the known couplings existing between the dopant and point defects (Self interstitials and Vacancies). All the charge states experimentally established for both the point defects and dopant/defects pairs are considered. The local equilibrium is not assumed for the various formation/dissociation reactions, thus the full system of continuity equations for both dopant, point defect and clusters is dynamically solved. Moreover, a model has been implemented for dopant precipitation when the solid solubility limit is reached.

The second component of the model deals with defect clustering effects where some of residual excess point defects nucleate and evolve into various forms of interstitial clusters (IC) like {311} defects, dislocation loops or vacancy clusters (VC). This evolution is explained by a competitive growth and dissolution of interstitial or vacancy clusters. These phenomena are described by the Ostwald ripening theory based on the reduction of free energy per defect of the extended defect which leads to the growth in size of the clusters during annealing [8] [9]. This phenomenon drives the interstitial super saturation and therefore modifies the dopant diffusion and clustering behavior.

Eventually, another part of excess point defects will be trapped by dopant atoms, to form immobile complexes, which makes up the third part of the model. In case of boron implantation, the model takes into account the formation of Boron Interstitial Clusters (BIC), like BI2+, B2I0, B3I2 + B4I3 +, where charge states are also considered according to [10]. While in the case of arsenic implantation, As atoms aggregate preferentiall with vacancies to form As2V0 or As4V0 complexes. This clustering phenomenon leads to an unwanted inactivation and immobilization of the implanted species.

Moreover, this simulation tool is flexible and expandable since other phenomena which may affect the dopant diffusion, i.e. the presence of carbon or fluorine atoms, can be easily incorporated by simply adding new chemical reaction and its coefficients to the model file. The use of such a complex physically based model involves many parameters which can be accurately determined only by a calibration methodology presented hereafter.


Calibration Methodology

The diffusion model presented above has up to 40 physical parameters including formation energies, kinetic coefficients, etc. Having in mind to keep them physically significant, an original calibration methodology has been developed [11]. It is based on specific statistical strategies, which include screening and sorting of relevant parameters, Design of Experiments (DoE) and Responses Surface Models (RSM) optimizations. This has been applied sequentially to all three components of the diffusion model. This calibration methodology has been successfully used on a broad range of experiments. Relevant simulations are presented below.


Simulation Results

A complete set of simulations has been performed using both BCA implant and advanced diffusion models. First, the interstitial clusters part of the model has been validated using the experimental data obtained by Cowern et al. [12] (Figures 1 and 2). Briefly, this experiment consists in observing the diffusivity of two boron marker layers after a silicon implantation at 40 keV to a dose of 2.1013at.cm-2. As exhibited in figure 1b, the model accurately predicts the two stages of the diffusion acceleration: the first plateau appearing at the earliest stage of the annealing, governed by the formation of small interstitial clusters, the second one characteristic of the competitive growth of {311} defects. Secondly, the calibration procedure was applied to tuning of the Pelaz experiment [13] in which a buried boron layer with a peak concentration above the 1019cm-3 level is annealed at 800oC for 35 min in the presence of defects created by a 2.1013cm-2 Silicon implant. It can be seen in Figure 3 that substantial portion of boron concentration may be kept immobile and inactive even below the solid solubility level, because boron atoms are trapped in mixed clusters (Boron Interstitial Clusters). The two previous experiments were used to calibrate the diffusion model, while its verification was carried out by characteristic implant/diffusion steps used in advanced technologies.

Figure 1. Simulation (solid line) of Cowern experiment (symbols) [12] where two boron MBE-pics diffuse after a Si implantation at 40keV/ 2.1013 cm-3. During the annealing, interstitial clusters will evaluate following through an Ostwald ripening scheme and will control the acceleration of the dopant diffusion.


Figure 2. Comparison between experimental (symbols) and calculated (solid line) I supersaturation evolution in Cowern experiment [12]. The good agreement between experiment data and simulations prove that the new model is able to handle all defects interactions and to simulate the various regimes in the dopant diffusion acceleration.


Figure 3. Comparison between experimental (symbols) and simulated (solid line) B profiles evidencing the influence of BIC’s. Data are extrated from Pelaz et al. (1997) [13].



As a validation test of the ability of our diffusion model to handle all the interactions existing during the dopant diffusion, three most representative implantation/diffusion experiments were chosen for simulations. The first experiment presented in Figure 4 shows that in the case of medium dose 20 keV boron implant annealed for 30 min at 800oC [14] most of the dopant concentration stays immobile because of the BICs formation. Moreover, the model provides a very important technological insight about the dopant inactivation, since a substantial portion of boron remains inactive. The last two examples presented in figures 5 and 6, simulate very low energy implant steps at 2 keV, respectively for boron [15] and arsenic [16], followed by a typical modern rapid thermal annealing (RTA) step. Complexity of these two simulations lies in the steepness of the implanted dopant concentration and the vicinity of the surface, where recombination phenomenon occurs. Unlike the boron experiment, where most of the dopant is active after the annealing, in the case of the arsenic the solid solubility model seriously affects both the shape of diffused profile and its inactive portion.


Figure 4. Comparison between experimental (symbols) [14] and simulated (solid line) B profiles implanted at 20 keV/2.1014cm-2, annealed at 800oC for 30 min. At high concentration, boron stay inactive due to the formation of Boron Interstitial Clusters (BICs).


The two last simulations (Figures 5 and 6) have been performed in 2D then compared to experimental data [15][16]. Finally a complete NMOS transistor has been fully simulated in 2D (Figure 7). During the annealing step, evolution and distribution of the <311> defect can be monitored (Figure 8) in order to better understand and thus optimize the dopant diffusion acceleration by defect engineering.


Figure 5. Simulation of ultra low energy implantation of B followed by RTA. (see ref. [15]). (Experimental data are represented with various symbols; lines are the simulation results).


Figure 6. Simulation of an ultra low energy implantation of As followed by spike annealing. (see ref. [16]). (Experimental data are represented with various symbols; lines are the simulation results).


Figure 7. Simulation of boron 2D distribution implanted at 2 keV/1014cm-2 and annealed at 950oC for 10 sec. Comparison between experimental data (symbols) (extracted from [15]) and profile extracted from the 2D simulation (solid line) evidences the quality of the new diffusion model.


Figure 8. 2D simulation of arsenic diffusion implanted at 2 keV/5.1014cm-2 and spike annealed at 1000oC with a ramp up estimated at 100oC/s. Most of the arsenic precipitate and stay inactive above solid solubility limit. (Experimental data are extracted from [16]).


All the coupled equations of the model are solved in a new specific solver called Difsim. One of the particularities of this solver is that it is able to integrate quickly new equations just by adding new reactions in the model files. Since this solver has been particularly optimized for this model, simulation times are short. For instance, it has only taken less than 7 min on a 2GHz-LinuxPC to simulate the 2D simulation of Figure 9.

Figure 9. 2D simulation of P-MOS transistor with boron diffusion after a low energy implant and RTA and with the presence of arsenic halo.


Figure 10. Evolution in time of the <311> defects concentration where at the early stage of the annealing they grow by exchange their interstitials to form clusters and finally dissolve according to the Ostwald ripening theory.



All the presented physical phenomena are successfully simulated with the BCA implantation module and with the new advanced diffusion model implemented within the SILVACO’s framework. A global and original calibration methodology has been used to provide the model with a physical set of parameters which allows to successfully simulate a very broad range of experimental setups, from low to high energy and dose implants followed by anneals ranging from a few seconds to tens of minutes. These results demonstrate that the simulation needs for key process steps of current technologies are fulfilled.


This article was presented at the E-MRS Spring 2005 conference. Full proceeding available at the Elsevier site: www.sciencedirect.com.

[1] Wang, et al., Phys. Rev. A44, p.1768, (1991)
[2] Azziz, et al., Physica Status Solidi 142, p. 35, (1987).
[3] Dort, et al., Solid-State Electronics 37, p.411, (1994).
[4] R. Smith Editor, Atomic and Ion Collisions in Solids and at Surfaces, Cambridge Univ. Press, (1997).
[5] Snyder and Neufeld, Phys. Rev. 97, p.1636, (1955).
[6] F. Boucard Ph.D. thesis (2003).
[7] Mathiot et al., J. Appl. Phys. 55, p. 3518.(1984).
[8] Claverie et al. Nucl. Instr. Meth. B 147, p.1 (1999).
[9] Ortiz et al., MRS Symp 669, p 5.6 (2001).
[10] Lenosky et al., Appl. Phys. Lett., 77, p.1834, (2000).
[11] F. Roger Ph.D. thesis (2002).
[12] Cowern et al. Phys. Rev. Lett. 82, p. 4460, (1999).
[13] Pelaz et al., Appl. Phys. Lett., 70, p. 2285, (1997).
[14] Solmi et al. J. Appl. Phys. , 88, p. 2135, (1991).
[15] B. Colombeau Ph. D. Thesis (2001).
[16] P. Fastenko Ph. D. Thesis (2002).
Download pdf version of this article.