# 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

*are:*

**ATHENA**- The spatial discretization scheme employed in
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.**SSuprem4** - The existing diffusion code may run into convergence problems during simulations on structures that involve polysilicon interfaces.
- The methods for linear system solving in
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 .**SSuprem4** - 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.

Dopants

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

*except in situations where the diffusion modeling of the latter is either wrong or artificial. Such exceptions include*

**SSuprem4**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**SSuprem4's**- A correction of the implementation of the fermi diffusion
model for indium, which was flawed in
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.**SSuprem4**

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,

*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.*

**ATHENA**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 |