**A Novel Approach to Three-Dimensional Semiconductor Process Simulation:
Application
to Thermal Oxidation**

*VASILY SUVOROV, ANDREAS HÖSSINGER,
ZORAN DJURIC and
NEBOYSHA LJEPOJEVIC
Silvaco Technology Center, Compass Point, St Ives, Cambridgeshire, PE27 5JL,
UK.
Journal of Computational Electronics Volume 5, Number 4, 2006c
2006 Springer
Science Business Media, Inc. Manufactured in The Netherlands*

**Abstract**

The paper presents a new approach to three-dimensional semiconductor process simulation that overcomes the problem of moving boundaries and mesh generation. Contrary to using unstructured meshes, the approach makes use of the level set method on fixed Cartesian meshes. A concept of multi-layer structure is introduced to capture an arbitrary complex structure. To handle a big geometrical scale ratio in a structure, the concept of adaptive mesh refinement is used. A special in-house finite-difference scheme is designed to approximate the relevant equations near material interfaces. In the bulk of regular nodes the standard finite difference schemes are used. Application of the approach to the modeling of oxidation of some typical types of structures used in semiconductor technology is demonstrated.

**Keywords**

Semiconductor technology, modeling, thermal oxidation of silicon.

**1. Introduction**

The 2005 edition of International Technology Roadmap for Semiconductors estimates that the technology-development costs can be reduced by 40% by 2013 due to use of Technology CAD. To achieve this goal, one needs to overcome a problem of mesh generation that has become a major issue in simulation of 3D process, mainly because of moving material interfaces. The standard approach that uses unstructured meshes involves very complex and often unstable mesh generation algorithms. Alternatively, we apply the level set approach [1]. All interfaces are given implicitly as signed distance (level set) functions, defined on fixed Cartesian meshes. Any changes of interfaces, caused by a given process, can be calculated by updating their level set functions. Every signed distance function satisfies the level set equation, which can be solved by standard methods.

The Cartesian meshes for the level set functions also serve as a framework for solving all other differential equations, which describe our simulated processes. The advantages of such an approach are numerous: numerical schemes for various equations are well known, stable and accurate; the coding is much simpler; the memory requirements per computational cell are smaller than for the finite elements methods. However, there are two problems associated with this methodology: approximating of general-type boundary conditions near interfaces and handling of a big scale ratio of sizes in complex structures. In recent decades, considerable progress has been achieved in resolving both issues (e.g. [2], [3], [4] ). In our development we use our original in-house numerical method, the method described in [3] and the concept of Adaptive Mesh Refinement (AMR) [4]. The above principles are cornerstones of a software package Victory-Process, Silvaco’s three-dimensional process simulation framework. In this work we demonstrate Victory’s capability to simulate numerically the most demanding process, namely thermal oxidation.

**2. Modeling of Thermal Oxidation**

Three different processes occur simultaneously during thermal oxidation: i) diffusion of oxidant through the oxide; ii) a chemical reaction at oxide/silicon interface which consumes silicon and generates volume expansion; iii) the deformation (flow) of the entire structure according to the rheological behavior of each material. The aim of the modeling and simulation is to explain and predict the resulting shape of the oxide and mechanical stresses developed in a structure during the process.

In order to represent a structure containing arbitrary number
of different materials and having arbitrary geometry, we developed the concept
of *multi-layer
structure*. Multi-layer structure is a set of signed distance functions
(*x; y; z*) defined in all computational domain. The level (*x;
y; z*) = 0 describes
an interface between a material (that is below the interface and < 0)
and the structure that is above the material (> 0). Moreover, the interface
spans the “vertical walls” of computational domain (i.e. the planes
x = 0 and x = a, y = 0 and y = b). As a result a material that is inside the
structure is normally outlined by two interfaces.

An elementary oxidation time step t is decomposed into a sequence of four simulation steps (2.1 - 2.4).

**2.1
Oxidant diffusion through the oxide**

At *time* = t, given all the interfaces and
the ambient oxidant flux, the distribution of the oxidant concentration *C*(*x,
y, z*) in the oxide has to be obtained by
solving the diffusion equation

(1) (2)

Where *D* is the diffusion coefficient, *k* is
the reaction rate, h is the gas-phase mass-transfer coefficient,* C*_{0} is the
equilibrium bulk concentration
in the
oxide and n is the normal to the correspondent interface.

Since our approach uses Cartesian meshes, we have developed a special finite-difference numerical scheme to approximate the diffusion equation at the points where material interfaces cross the regular Cartesian grid. In order to apply the boundary conditions (2) at such ‘irregular’ points, we are adding cross-points to the problem formulation. Whenever a mesh line of the Cartesian mesh crosses an interface, a crosspoint is introduced as illustrated in Figure 1. Once these cross-points are inserted, one of three types of equations can be formulated for each Cartesian node (including cross-points). Firstly, the standard finite-difference scheme can be used for all points inside the oxide. Secondly, nearly standard finite difference equations can be applied for all regular nodes which are close to the interface. The only difference is, that for those equations the cross-points have to be taken into account instead of neighboring regular nodes. Thirdly, only the cross points need a special treatment. While it is straightforward to apply Dirichlet boundary conditions, a special interface equation has to be derived for the case of Neumann or mixed boundary conditions, which takes into account tangential derivatives. One derivative of ∂C/∂n is always directly accessible, since each cross-point is sitting on a mesh line. The other two derivatives can be obtained by also considering points in tangential direction. Tangential direction means the direction from the cross-point to a neighboring cross-point.

Figure 1: Discretization on the Cartesian mesh near arbitrary interfaces. |

Additionally there is one special case when the interface is passing very close by a regular node or exactly through a regular node. In this situation, the standard discretization for ∂C/∂n is used which takes only regular nodes into account. Thereby, numerical instabilities, due to extremely large coupling coefficients, can be avoided.

**2.2 Propagation of the oxide/silicon
interface**

After solving the diffusion equation the time is set to *time* = t
+ t. While the time dt elapses, a thin layer of silicon will be burnt. To
find a new
position of the oxide/silicon interface the level set equation is solved
for the distance
functions that represent the interface oxide/silicon.

(3)

where N is the number of oxidant molecules incorporated into
unit volume of the oxide, is *Si* *SiO*_{2}
expansion coefficient. Equation (3) is solved by methods which use the high-order
nonoscillatory schemes for
Hamilton-Jacobi
equations [5]. The space derivatives are approximated by essentially nonoscillatory
schemes up to fourth order; all terms which depend on space derivatives are
approximated by monotone fluxes of various types: Lax-Friedrich, Godunov,
Roe
with entropy correction. As expected, the least diffusive, and best suited
for our simulations, turned out to be Godunov numerical flux. For updating
equation (3) in time, a class of TVD Runge-Kutta discretization was used,
and the implementation includes up to third order schemes. The entire numerical
procedure is explicit, and therefore the CFL condition puts the limitation
on the maximum time step to be used. However, the numerical solver is very
efficient and accurate.

**2.3 Calculation of the volume expansion resulting
from the chemical reaction**

The layer of the burnt silicon turns into a thicker
layer of oxide with a significantly larger volume that forces all the overlaying
structure to deform.
To find the
velocity field which describes this deformation we solve the creep equation

(4)

where *S _{ij} * is a Cauchy stress tensor. For
incompressible fluid
the tensor can be decomposed into two parts:

(5)

where * _{ij}* is the shear stress tensor defined by the rheological
model of a material. The commonly assumed model for the oxide and nitride is
the Maxwell
viscoelastic fluid with the following constitutive equation

(6)

where *u _{i}* is a velocity component. The relaxation
time = µ/

*G*gives a reference time that separates two different flow regimes: for t << the liquid behaves similar to linear elastic body, for t >> the liquid exhibits viscous behavior. For the examples shown in this paper we use the viscous model both for the oxide and nitride, neglecting the second term in the left-hand side of the equation (6). The further simplified assumption is that the oxide and nitride viscosities are stress-independent. The model can be justified for high temperatures (T 1000

^{o}C) when µ and are relatively small. In this case the system (4-6) becomes the system of the Stokes equations with the following boundary conditions:

(7)

(8)

Here [·] denotes the jump across any liquid/liquid interface; is the surface tension coefficient and is the surface curvature. We note that the oxygen atmosphere can formally be regarded as an incompressible liquid with µ = 0.

Dirichlet boundary conditions (7) are implemented by extrapolating the given velocities values from the cross-points to the regular mesh lines. In this case one can use the standard finite-difference discretisation schemes.

Another problem we had to overcome when solving the system (5-8) is the problem of discontinuous viscosity: viscosity has a jump across oxide/nitride or oxide/oxygen interfaces. In order to resolve the problem we used the method suggested in the paper [3], where viscosity near the interfaces is smeared out to the underlying finite difference grid using regularized Heviside function.

**2.4
Propagation of the interfaces in the deformation velocity field**

The final step
updates all the liquid/liquid interfaces that are above the silicon/oxide one
by solving the advection equation for all correspondent
level set functions in the multi-layer structure representation.

The same numerical schemes are applied for this numerical step as for the interface propagation modeling the conversion of silicon into silicon dioxide.

**3. Results**

For illustration purposes we have modeled three typical
oxidation process steps: trench oxidation, LOCOS, and polysilicon oxidation.
Figures 2, 5 and
8 show
initial structures produced by Victory’s etching/ deposition module.
Figures 6 and 9 demonstrate two different types of three-dimensional bird’s
beak effects. Pictures 3, 6 and 9 presents the familiar ‘2D’ pressure
contours at the edges. As expected, the pressure is high near curved surfaces:
it has positive values (red, compressive) near a concave surface and has
negative values (blue, tensile) for a convex one. The pressure is measured
in the values
of *P _{ref} *= 7.5 · 10

^{8}Pa. Three-dimensional isosurfaces of the pressure values that are close to the maximum ones inside the structures are presented in the Figures 4, 7 and 10.

Figure 2: Input structure for the trench oxidation generated by the etching/deposition module. |

Figure 3: The trench structure and the pressure distribution after 4.4 minutes of oxidation. |

Figure 4: Iso-surfaces of tensile pressure P = -2070 Mpa (blue) and the compressive one P = 5250 MPa (red) in oxide. Yellow surface is the oxide/silicon interface. 4 |

Figure 5: Input structure for the locos process built by Victory’s etching/deposition module. |

Figure 6: The locos structure and the pressure distribution after 4.3 minutes of oxidation. |

Figure 7: Iso-surfaces of the tensile pressure P = -1125 Mpa (blue) and the compressive one P = 975 MPa (red) in oxide. Yellow surface is the oxide/silicon interface. |

Figure 8: Input structure for the reoxidation of two oxidizable regions |

Figure 9: The polyreoxidation structure and the pressure distribution after 5 minutes of oxidation. |

Figure 10: Iso-surfaces of the tensile pressure P = -4500 MPa (blue) and the compressive one P = 27000 MPa (red) in oxide. Yellow surface is the oxide/silicon interface. |

**4. Conclusion**

The state-of-the-art finite-difference methods for the simulation of thermal oxidation is presented which can handle arbitrary complex geometries. This enables us to design a new generation process modeling tool that avoids meshing/re-meshing procedures.

**References**

- J. A. Sethian, Level set methods and fast marching methods, Cambridge University Press (1999).
- S. Deng, K. Ito, and Z. Li, Three-dimensional elliptic solvers for interface problems and applications, J. Comp. Phys. 184, 215 (2003).
- M. Sussman, P. Smereka, and S. Osher, A level set approach for computing solutions to incompressible two-phase flow, J.Comp. Phys. 114, 146 (1994).
- P. Colella et al. Chombo software package for AMR applications: design document, available at http://seesar.lbl.gov/anag/chombo/ .
- S. Osher and C. Shu, High order essentially nonoscillatory
schemes for Hamilton-Jacobi equations, SIAM J. Nubmer.
Anal. 28, 907 (1991).