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

Fast Diffusion Module

A new diffusion algorithm based on a Galerkin method with linear finite elements, an extremely fast sparse matrix solver, and object oriented physical modeling is one of the new features implemented in ATHENA version 5.0. This module is an alternative to the existing code, thus providing the user with a choice between using the old diffusion module and the new algorithms.

The Implicit Linear Finite Element Method (ILFEM) module is available for diffusion and oxidation steps. Switching between the two modules is done using a method statement in the input file, giving the user a flexible choice between numerics for each individual diffusion run.

The advantages in adding a new diffusion module to ATHENA are:

  • The spatial discretization scheme employed in SSuprem4 causes a slow-down during simulations on grids with poor geometry. Due to the inability of the first order finite volume method to deal with triangles with large aspect ratios. An excessive amount of time is spent doing Lawson flipping, whenever skinny triangles start to accumulate at the moving Si-SiO2 interface during oxidation.
  • The existing diffusion code may run into convergence problems during simulations on structures that involve polysilicon interfaces.
  • The methods for linear system solving in SSuprem4 were inherited from the original Stanford SUPREM-IV[1]. There are not fully up to date with most recent progress within the field. For optimal efficiency, most contemporary simulators use iterative sparse linear solvers, whenever possible1 .
  • When a high concentration of <311> interstitial clusters are present in the structure there may be numerical problems.
  • Improvement modeling of interstitial traps, as described by P. Griffin [2].
  • The code on which the original diffusion simulator is based, was created before object-oriented programming paradigms and tools to support it became prevalent. When programs based on procedural programming grow large, they tend to become unmanageable, thus degrading both maintainability as well as extendibility. For rapid implementation/ testing of new physical models, well organized base code is of paramount importance.

In the following sections, the new module will be described in more detail, especially how it addresses some of the above mentioned problems in SSuprem4.

The Physics

It may be shown that it is not energetically favorable for an impurity atom to migrate on its own in a semiconductor lattice. Overcoming the large Coulomb repulsion in the diffusion saddle point makes this mechanism rather prohibitive, so dopant diffusion is believed to predominantly take place through interaction with point defects[3]. A dopant atom can form a pair with either an empty lattice site (a vacancy, V) or a substrate atom sitting between regular lattice sites (a self interstitial, I), and the pair as a whole then diffuses through a number of inversion cycles. The set of reaction-diffusion equations solved in ATHENA is based on such a dopant/point-defect pair diffusion mechanism[1]. In the following we ignore the fact that dopant-defect pairs can exist in different charge states.


The continuity equation for dopants is where CA designates the total concentration of dopant A
and JAX is the pair flux. Notice that this expression explicitly states that dopant atoms can only move via the pair mechanism. The flux expression for dopants is

where DAX is the pair diffusivity, CA+ is the electrically active dopant concentration, C*X designates the equilibrium defect concentration, and Z is the dopant charge number (+1 for donors and -1 for acceptors). The dopant flux expression (2) incorporates contributions from both the dopant gradient, a supersaturation of point defects, as well as an internal electric field. The pair diffusivity is actually a composite entity, taking into account the distribution of pairs over different charge states as well as fermi level enhancement effects, according to the following expression

where c is the charge state, n and ni are the actual and intrinsic carrier concentration, respectively, and the temperature dependency of the intrinsic pair diffusivity, is given by an Arrhenius expression,

In the fermi model, vacancies and self interstitials are assumed to be in thermodynamical equilibrium, and the defect supersaturation ratio becomes unity in (2).

Point Defects

Due to implantation damage and injection during oxidation and nitridation, the point defect populations can in general not be assumed to be in equilibrium, for which reason explicit representation and evolution of these species is needed.The continuity equations for vacancies and interstitials are

where R is a recombination term that accounts for the creation/annihilation of interstitial-vacancy pairs in the semiconductor bulk,

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

Diffusion Modeling

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.


Example Model Nodes ATHENA 4.3 ILFEM Speedup
MOSFET two.dim 2829 583s 240s 2.43
Vertical DMOS fermi 458 119s 99s 1.2
Power Device fermi 8250 318s 85s 3.74
Arsenic Activation full.cpl 672 429s 68s 6.31