On continuum simulations of the evolution of faulted and perfect dislocation loops in silicon during post-implantation annealing

Predictive models of the damage evolution during post-implantation annealing are important for predicting the performance of final devices. An existing model of the thermal evolution of dislocation loops in silicon during post-implantation annealing has been recalibrated to better capture the distinction between faulted and perfect dislocation loops. The calibration was based on experimental findings from the literature on the simultaneous presence of both faulted and perfect dislocation loops after post-implantation annealing.


Introduction
During post-implantation annealing, silicon self-interstitials (I) and silicon vacancies (V) generated during the implantation process are assumed to recombine, and the excess of Is is assumed to start forming clusters. For certain conditions, extended defects such as {311}-defects and dislocation loops, faulted (FDL) and perfect (PDL), are formed. This is described by, e.g., Claverie et al. [1]. The dominating types of defects formed depend on the implantation damage and the annealing temperature and time. The damage evolution during annealing contributes to an enhanced dopant diffusion because of the resulting supersaturation of Is, and larger defects, when they form, may result in higher leakage currents in the final device. Dislocation loops are associated with higher implant energies and doses, and higher thermal budgets, used in the manufacturing of, e.g., power electronic devices and solar cells.
Computationally efficient simulation models, which describe the damage evolution, exist in the literature [2][3][4], to mention a few. The state-of-the-art model by Zechner et al. [2] describes the evolution of small I-clusters up to a certain size and uses a moment-based approach to describe {311}-defects. It was extended by Zographos et al. [3] by including a moment-based method for dislocation loops (without distinguishing between FDLs and PDLs). Wolf et al. [4] found that the model by Zographos et al. [3] resulted in too high dissolution velocities of loops at high temperatures and a diverging mean loop radius for long annealing times at high temperatures, contradicting experimental findings [5]. To account for that, they extended the model by Zographos et al. [3] to include PDLs and a different model approach for the evolution of the mean loop radii, implementing the saturation of Ostwald ripening [5]. For its calibration, presented in Wolf's thesis [4], the focus was predominantly on conditions relevant for solar cell processing (long high-temperature anneals), where they mainly considered the sum of dislocation loops rather than the distinction between the two types.
In this work, the focus was to recalibrate the model by Wolf et al. [5] to better describe measured data where FDLs and PDLs coexist after post-implantation annealing. The better the individual defect types can be described, the clearer the overall picture will be, contributing to better predictions of enhanced dopant diffusion and performance deterioration due to leakage currents in final devices.

Simulation model
The model of Wolf et al. [5] is presented here to give an overview of the available tuning parameters. The thermal evolution of dislocation loops of type L is modeled by 1 3 It is expressed analogously for FDLs and PDLs. C FDL , C PDL , and C 311 denotes concentrations of Is contained in the respective loop types and in {311}-defects, and n FDL , n PDL , and n 311 , the respective loop densities. k 311 L governs the {311}-defect-to-loop transformation together with the scaling factor k 311 L . For example, k 311 L = 0.5 has the meaning that {311}-defects twice the size of their mean size transform into loops [3]. k I L governs the loop growth from Is. C I and D I are the concentration and diffusion coefficient of Is. The mean loop radii are denoted r FDL and r PDL . FDL ∕r FDL and PDL ∕r PDL represent the normalized variances of the loop size distributions, which are assumed to be time independent [4]. C * I,L (r L ) , the concentration of Is in the vicinity of a loop of type L and size r L , is given by where k is Boltzmann's constant. ΔE L for FDLs and PDLs are given by and (1) (3) respectively. is the stacking fault energy, Ω is the atomic density of silicon, b FDL and b FDL are the burgers vectors b 2 PDL = 3 2 b 2 FDL , G is the shear modulus. A L ∈ [−2.77, 1] represent the kernel energies of the respective loops [5].
The mean loop radius for loops of type L is estimated according to assuming circular loops (discs) [5], where d 111 is the atomic density of {111} planes ( 1.57 × 10 15 cm −2 ).
The complete system of equations also includes the mentioned model by Zechner et al. [2], bulk recombination of I and V, and boundary conditions, already implemented in Sentaurus Process of Synopsys [6], the software used for all simulations presented in this work.

Simulation results
In this section, the reference data are presented together with the corresponding simulations results of the implantation and the post-implantation annealing. Simulation results using the parameter set of Wolf [4] are shown to motivate the need of a recalibration, and finally, the recalibration is presented.

Experimental details
The experimental data used as reference for the calibration, taken from Cristiano et al. [7] and Pan et al. [8], are shown in Fig. 1. Each experiment is represented with two plots, one for the total number of loops ( n tot L ) vs. time (top) and one for the mean loop radii assuming circular discs vs. time (bottom). The experimental details are summarized in Table 1. In Ge1 and Ge2, Ge was implanted with a dose of 2 × 10 15 cm −2 and an energy of 150 keV, resulting in 175 nm thick amorphous layers [7]. In Ge2, the amorphous layer was reduced to 45 nm through etching by anodic oxidation. Both samples were annealed by rapid thermal processing (RTP) for 10-400 s in N 2 . In Si1 and Si2, Si was implanted with a dose of 2 × 10 15 cm −2 and an energy of 50 keV, resulting in 100 nm thick amorphous layers [7]. Both samples were first annealed for 10 min at 900 • C in N 2 in a conventional furnace, followed by annealing at 900 • C for an additional 15-120 min in N 2 for Si1, and in O 2 for Si2. In Si3 and Si4, Si was implanted with a dose of 1 × 10 16 cm −2 and an energy of 50 keV [8]. The thicknesses of the resulting amorphous layers were not reported in the paper. Si3 was annealed in a conventional furnace at 850 • C for 5-960 min in N 2 . Si4 was annealed by RTP at 1000 • C for 10-400 s in N 2 .
Unfortunately, the ramp-up rates of the considered experiments were not stated, only whether a conventional furnace or RTP was used. For the furnace anneals, a ramp-up rate of 10 K/min was assumed, starting at 500 • C . For the anneals with RTP, a ramp-up time of 5 s from 500 to 900 • C was assumed, followed by a ramp-up time of 3 s to reach 1000 • C in the Si4 case. This is based on the ramp-up times specified in the work of Bonafos et al. [9] (from the same research group as Cristiano et al. [7]).

Implantation damage
Before simulating the post-implantation anneal, it is crucial to have a good description of the implantation damage. In this work, the implantation step was simulated with the Monte Carlo simulator of Sentaurus Process [6]. First, the effective +n model based on the work of Hobler and Moroz [10] was used to describe the implanted damage. The +n represents the number of remaining Is per implanted ion beyond the a/c interface. However, the calculated total number of Is in loops based on the measured data after annealing in inert ambient is higher than the simulation results directly after implantation. This motivates the use of a larger +n than what was calculated by the effective +n model.
In the simulations, the thickness of the amorphous layer is determined from the concentration of displaced atoms. If the concentration exceeds the amorphization threshold, the corresponding region is assumed to be amorphous. For the experiments of Cristiano et al. [7], the simulations predicted amorphous layers of approximately 175 nm and 100 nm for the Ge and Si implantation, respectively. For the experiments of Pan et al. [8], the simulation results predicted a ca. 130 nm thick amorphous layer. The default amorphization threshold from Sentaurus Process [6] was used for the simulations. The simulation results for the considered experiments are summarized in Table 2, reporting also the used +n for each implantation.

Simulations with Wolf's parameters
To reproduce and verify the simulation results of Wolf [4], it is necessary to also have the same starting conditions, not only the same model parameters for the damage evolution. Wolf [4] described the simulation results of the implantation with the a/c depth and the net number of Is after implantation. Based on that, his results were reproduced by adjusting +n and the threshold for amorphization. For Ge1 and Ge2, the listed information in Wolf's work [4] for the Ge-implantation experiments of Bonafos et al. [9] (with the same specification) was used.
The simulation results of the post-implantation annealing are shown in Fig. 1 as dashed lines. In Fig. 1a-b, the PDLs dissolve according to the experimental data, especially in the Ge2-case. In the simulation results, n tot L and the corresponding r L saturate for both loop types, and r PDL >r PDL . In Fig. 1c-d, n tot FDL and n tot PDL saturate, like the experimental data indicates. Although, n tot PDL is significantly underestimated for both Si1 and Si2. For Si2, r FDL is in good agreement with the data, but the decrease in r PDL at ca. 80 min is not reproduced. In Fig. 1e, for Si3, n tot PDL is underestimated by orders of magnitude, and r FDL and r PDL are both overestimated. In Fig. 1e, for Si4, n tot L is overestimated for both loop types. r FDL is in good agreement with the data, but r PDL is underestimated.
The formation energy of k 311 PDL from Wolf [4] was adjusted from 6 eV to 6.50 eV for these simulations because assuming 6.0 eV from the precision given in the original work resulted in a poor agreement with the therein reported simulation results of Si1 and Si4 (the remaining experiments were not considered by Wolf [4]). It should be mentioned that there is a slight mismatch between original work and the simulation results shown here. Whether the difference is due to the precision of the listed parameters, different assumptions about the implanted damage or the ramp-up rates, the implementation, or a combination of these, is hard to say.

New calibration
The calibration task is challenging. For example, it is hard to describe Ge1 and Si1 with one parameter set. In Ge1, the PDLs dissolve within ca. 100 s, but not in Si1, even though the annealing temperature was the same. Cristiano et al. [7] motivated the different findings with the difference in implantation damage. In the Si1-case, there is a higher excess of Is in the end-of-range region after implantation compared to the Ge-cases. Another difference is the ramping rate, the ramp-up phase is much longer for the relatively slow ramping of a conventional furnace.  Table 1: a Ge1, b Ge2, c Si1, d Si2, e Si3, and f Si4. The data in a-d and in e-f were taken from Cristiano et al. [7] and Pan et al. [8], respectively The parameters were calibrated by careful hand-tuning. The focus was to increase n tot PDL at the lower temperatures, while maintaining a good description of the radii evolution. To evaluate the error, a loss function including weighted errors in the number of loops (on a logarithmic scale) and the radii of the loops was used. The parameter values are summarized in Table 3 and the corresponding simulation results are shown as solid lines in Fig. 1.
A FDL was chosen as 0, just like in previous works [3][4][5] and A PDL was set to 0 as well, compared to − 2.65 [4]. This basically means that the PDLs are assumed to be less stable compared to before. The new value is in line with the arguments presented by Cristiano et al. [7], regarding the stability inversion, which was not maintained by A PDL = −2.65.
Following the arguments of Wolf [4], C * I,L (r L − L )− C * I,L (r L + L ) should be larger for PDLs than for FDLs for all loop sizes, representing the more pronounced Ostwald ripening for PDLs [8]. This is achieved with the new values for PDL ∕r PDL and FDL ∕r FDL . Just like in Wolf's work [4], PDL ∕r PDL > FDL ∕r FDL , but both new values were chosen slightly smaller.
The reaction rate constants k 311 FDL and k 311 PDL were both chosen larger compared to Wolf [4], with k 311 FDL < k 311 PDL , resulting in n tot PDL closer to the measured values. k 311 FDL,PDL were not changed based on the argumentation of Wolf [4]. Finally, the new values for k I FDL and k I PDL are significantly smaller at low temperatures compared to the previously suggested values. This change results in a slower growth of the mean loop radii at lower temperatures, especially significant in the ramping-up phase of the anneal.
The simulation results after calibration better describe the experimental data than the simulation results based on Wolf's work [4]. Although, in Fig. 1a-b, the dissolution of PDLs is still not captured and r PDL is still larger than r FDL after calibration. For Ge1, n tot FDL is in good agreement with the data, and r FDL is slightly underestimated. For Ge2, both n tot FDL and r FDL are somewhat underestimated. On a more positive note, the experimental sum of loops, in both Ge1 and Ge2, are well described by the simulation results.
In Fig. 1c, the simulation results much better predict n tot L for both loop types. The results also show a good fit for r FDL and r PDL . Similarly, in Fig. 1d, n tot L is in good agreement with the data for both loop types, but the decay in r PDL at longer annealing times is not captured.
In Fig. 1e, n tot L is slightly underestimated at longer annealing times for both loop types, and r FDL and r PDL are both overestimated. Nevertheless, the agreement is improved compared to the dashed lines. Finally, in Fig. 1f, n tot FDL , n tot PDL , and r PDL are well described by the simulation results, but r FDL is slightly overestimated.
In principle, the model can reproduce the Ge-implantation and Si-implantation experiments, but not with a constant parameter set. With the small number of relevant studies available, it is hard to analyze to what extent the experimental conditions are truly comparable. Another aspect that could contribute to the discrepancies are the assumptions about generation of the intrinsic interstitials during ion implantation, since the initial damage before the anneal influences its evolution. In order to get better overall description and assessment of the data, more experiments would be needed.

Conclusions
It was possible to better describe the considered experimental data with both FDLs and PDLs, by calibrating the model for the evolution of dislocation loops by Wolf et al. [5]. Although the fit is better, there is still room for improvement. More detailed information about the implantation damage and the annealing conditions, including ramping rates, would be beneficial for a better assessment of their respective impacts. Nevertheless, the new calibration contributes to the continued development of I-cluster models that have been used quite successfully for around 20 years to describe transient diffusion in silicon and SiGe. These models can always be improved further, of course, e.g., by considering the contribution of doping atoms.

3
need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.