New, Fast Numerical Algorithm for Diffusion Modeling Implemented in ATHENA Version 5.0

with KR a rate constant. The point defect flux expressions are given by

which correctly accounts for the effect of an electric field on the charged portion of the defects (the defect equilibrium concentrations, C*V and C*I, are fermi level dependent).For the sake of brevity, neither boundary conditions nor interactions with interstitial traps, <311> clusters, nor dislocation loops have been mentioned.

The presence of the pair fluxes, JAV and JAI, in (5) and (6) represent the contribution from paired defects to the total time rate of change of the defect populations, and establishes a two-way coupling between the diffusion of dopants and point-defects. Therefore, the three-species model defined by eqs. (1), (5), and (6) is called the fully coupled diffusion model. For simplification, summation over all present dopants and charge states has been omitted. In reality, there should be a pair flux correction term from every combination of these degrees of freedom. Under circumstances where the pair fluxes become very small in comparison to the fluxes of free point defects, the fully coupled correction terms can be omitted altogether. This is equivalent to zero binding between dopants and defects, and is called the two-dimensional diffusion model.

The system of coupled partial differential equations defined by the the equations (1), (5) and (6) is both non-linear and stiff, so special attention is required in the numerical treatment. Particularly, the time discretization method used for integrating a stiff PDE system must be robust in order to avoid numerical instability2. The original diffusion module uses a composite trapezoidal second order backward differential rule (TR-BDF2) for the time integration, while ILFEM uses a backward euler one step rule. The TR-BDF2 [4] method has a smaller local truncation error (third order in time) than backward euler (second order in time), which, on the other hand, is easier to implement. Both methods, however, do satisfy the stability requirements needed to deal with stiff problems.

The segregation of dopant impurities across interfaces is modeled by means of a simple first order kinetic model

where C1 and C2 are the particle concentrations in the immediate vicinity of the interface, h12 is the transport velocity, and M12 is the segregation coefficient, defined as the ratio between solubilities of the dopant in the two materials sharing the interface. For the new diffusion module, h12 has been set artificially large, thus forcing the dopant to instant equilibration on the interface. This penalty method, of course, does not lend itself to the modeling of transient segregation phenomena, but, on the other hand, is an efficient method for the damping out of numerical oscillations that tend to originate on the interface where the dopant profile is discontinuous.


Spatial Discretization

The finite volume method, runs into difficulties when obtuse or skinny triangles are present in the grid. The method depends on the correct calculation of Voronoid cell areas as well as coupling coefficients, both of which fail in the partitioning of an obtuse triangle. This causes numerical havoc, inducing such unphysical effects as reversing the direction of diffusion flow through cells or distortion of the flux. To prevent the forming of triangles of bad aspect ratios during oxidizing runs, the grid topography has to be constantly monitored throughout the simulation, and when a bad triangle is encountered, a Lawson flip of edges is required. This extra overhead, not surprisingly, may cripple the performance during simulations of heat cycles in oxidizing ambients.

A standard Galerkin linear finite element methods has been chosen for the ILFEM diffusion module. This method is of second order precision in the spatial variables, uses linear test functions, and is well suited for irregularly shaped geometries. Most important of all, though, is the fact that it is totally immune to the pathological grid effects mentioned above.


Iterative Sparse Linear System Solver

The ILFEM diffusion module uses a very fast iterative sparse linear solver that employs domain decomposition preconditioning. A significant performance increase in linear solving has been observed in benchmark tests. The performance gain, though, depends on various factors such as the size of the grid, and the condition number and sparsity structure of the stiffness matrix. The relative performance is best for large, symmetrical matrices.


Object Oriented Physical Modeling

The new module is a total rewrite of the diffusion solver, using object-oriented design from the bottom-up. The objective has been to both reproduce SSuprem4's basic three diffusion models as well as to provide a solid frame work for future implementation of new physical models that lends itself to rapid development, robust code, and easy verification and validation.

The design incorporates the following features

• Models implemented as a C++ class hierarchy
• Dynamical creation of physical model objects (:created when needed)
• Smart model switching
• Dynamical creation of FEM matrix generators
• Localization of all model parameter information
• Logging of all used model parameters for debugging and control
• Extensive use of polymorphism
• Easy introduction of new materials and impurities through localization of control logic
• Separation of models, grid information, and numerics

From a modeling perspective the goal of the new diffusion module, as implemented in the 1999 ATHENA release, has been to reproduce SSuprem4 except in situations where the diffusion modeling of the latter is either wrong or artificial. Such exceptions include

• SSuprem4's unphysical pile-up of dopants at oxide/gas interfaces. Using the new diffusion module, the dopant profile will instead make a dip near the oxide/gas interface because of dopant loss to the ambient
• A correction of the implementation of the fermi diffusion model for indium, which was flawed in SSuprem4 is also included
• SSuprem4 modeling of the interstitial trap reaction, has much too high trap reaction rate. The rate has been changed to a more physically appropriate value, resulting in a more smooth fall-off of the edge of the interstitial profile during coupled diffusion.

For the current release four models have been ported to ILFEM: The fermi model, the two dimensional mode, the fully coupled model, and finally the high concentration fully coupled model as formulated by S. Crowder [5].

The following illustrations show 1D diffusion profiles for four commonly used dopants in different physical scenarios (dose, implant energy, diffusion model, anneal temperature). Each graph contains an initial profile and the results obtained through ATHENA 4.3 and ILFEM, respectively. All examples employ a simple SiO2/Si structure in order to model the diffusion of the dopant in silicon as well as its segregation across a silicon/oxide interface. For each case a 0.5 micron oxide layer was deposited after dopant implant, using 20 grid divisions. Notice that all profiles are displayed in a semi logarithmic plot.

Figure 1 shows a simulation of a one hour fermi diffusion at 900 C of a medium dose, high energy boron implant (5e14 cm-2, 80 keV). Since ni, the intrinsic carrier concentration is about 4e18 cm-3 at this temperature, the doping is extrinsic between 0.2 microns and 0.62 microns, for which reason some fermi level enhanced diffusion should be expected. The agreement between ATHENA 4.3 and ILFEM is seen to be overall good, ATHENA 4.3 being slightly more diffusive in the tail region of the profile. There also seems to be good consistency with respect to the amount of segregation of boron into the oxide layer. One should bear in mind, though, that ILFEM employs infinite high transport velocities across all interfaces, so during the time of transient segregation the profiles should be expected to differ slightly in the oxide.

Figure 1. Boron anneal using ILFEM and ATHENA.

Figure 2 shows a fifteen second fully coupled diffusion at 1000 C of a medium dose phosphorus implant (5e14 cm-2, 35 keV). The implant model used was dual pearson with a damage factor of 0.1. Since the intrinsic carrier concentration is about 1.0e19 cm-3 at 1000 C, the joint effects of fermi level enhancement and TED (Transient Enhanced Diffusion) are to be expected in the diffused phosphorus profile. The agreement between ILFEM and ATHENA 4.3 is seen to be excellent both below as well as above the characteristic phosphorus kink.

Figure 2. Fully coupled phosphorus diffusion using ILFEM.

Figure 3 shows the simulation of a one second 1100 C RTA profile of a high dose arsenic implant (5e15 cm-2, 80 keV), using the two.dim diffusion model. The intrinsic carrier concentration is about 1.2e19 cm-3 at 1100 C, so the silicon less than 0.3 microns from the SiO2/Si interface is highly extrinsic. Apart from a slight deviation in the peak region of the arsenic profile which might be due to a minor difference in the implementation of the AsV clustering model, the agreement between ATHENA 4.3 and ILFEM is again good.

Figure 3. Defect enhanced diffusion of arsenic using ILFEM enhanced.

Figure 4 shows a two hour fermi diffusion at 1100 C of a high dose antimony implant (5e15 cm-2, 80 keV). Again the interface region of the semiconductor is highly extrinsic, so fermi level enhancement is expected. Due to the low solid solubility of substitutional antimony in silicon, a peak of clustered/precipitated antimony can be observed near the silicon/oxide interface. While the agreement between ATHENA 4.3 and ILFEM is still good in the dopant peak region, there is a slight deviation in the tail region as well as in the oxide.

Figure 4. High temperature antimony diffusion in ILFEM.

Benchmarks Tests

Four different benchmark tests have been carried out in order to compare the performances of the two diffusion simulators. The first test involves an RTA activation of a high dose source-drain implant, the second test consists of the process part of the formation of a vertical DMOS device the third test is 420 minutes of fermi diffusion on a power device, and finally the fourth test simulates one second of fully coupled diffusion of a 1.0e15 cm-2 dose arsenic implant with high implant damage. All simulations were executed on a Sun UltraSPARC-IIi 300MHz workstation.

The results are presented in Table 1.