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₀ 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₂ 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^oC) 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⁸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

1. J. A. Sethian, Level set methods and fast marching methods, Cambridge University Press (1999).
2. S. Deng, K. Ito, and Z. Li, Three-dimensional elliptic solvers for interface problems and applications, J. Comp. Phys. 184, 215 (2003).
3. 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).
4. P. Colella et al. Chombo software package for AMR applications: design document, available at http://seesar.lbl.gov/anag/chombo/ .
5. S. Osher and C. Shu, High order essentially nonoscillatory schemes for Hamilton-Jacobi equations, SIAM J. Nubmer. Anal. 28, 907 (1991).
   

Download PDF of this article