Subsections

• 5.1 Simulation Flow of ELSA
• 5.2 Initialization
• 5.3 Transport of Particles
• 5.4 Extending the Speed Function
• 5.5 Fast Marching Method
• 5.6 Narrow Banding Technique
• 5.7 Combining Narrow Banding and Extension of the Speed Function
• 5.8 Coarsening Algorithm
• 5.9 Advancing the Level Set Function Using Finite Difference Schemes
• 5.10 Iterative Solver
• 5.11 Input File of Three-Dimensional ELSA
• 5.12 Fully Three-Dimensional Process Simulation
• 5.13 Stability of the Simulator
• 5.14 Summary

5. The General Purpose Topography Simulator ELSA

In this chapter the application of the level set and fast marching methods for the development and implementation of the topography simulator ELSA in two and three dimensions is presented. ELSA rest on many techniques, including a narrow band level set method, fast marching for solving the Eikonal equation and extension of the speed function, transport models, visibility determination, and an iterative equation solver.

The development and implementation of a two-dimensional topography simulator based on the level set method have begun at our institute by Dr. Heitzinger [35]. It has been written in Lisp. However, its use was only limited to the simulation of the deposition of layers into single features.

5.1 Simulation Flow of ELSA

Every feature-scale topography simulation must take the following three main steps into account. First, the transport of particles above the wafer in the boundary layer must be simulated. This can happen in the radiosity regime as described in Section 5.3. Second, the fluxes of particles at the wafer surface found in the previous step determine the chemical reactions that take place at the wafer surface. Third, the surface of the substrate moves according to the fluxes found in the second step.

Figure 5.1: Overview of simulation flow as combination of a physical transport model and surface evolution using the level set method.

The simulation flow is shown in Figure 5.1. The simulation stops when a prescribed time is reached or when a layer of prescribed thickness has been deposited. In the next sections the different techniques used during the implementation of ELSA are explained and discussed.

5.2 Initialization

In order to apply the level set method a suitable initial function $\varphi (0,\mathbf{x})$ has to be determined first. There are two requirements: first it goes without saying that its zero level set has to be the surface given by the application, and second it should essentially be a linear function so that linear interpolation can be applied in the final surface extraction step.

The signed distance function of a point from the given surface has shown to be a good choice for the initial level set function. This function is the common distance function multiplied by minus or plus one depending on which side of the surface a grid point lies. The common distance function of a point $x$ from a set $S$ is then defined by $d(x,S) := \inf_{y\in S} d(x,y)$, where $d$ is a metric, usually the Euclidean distance.

Consider a fixed orthogonal grid for a whole simulation domain. Because of the ease of two-dimensional illustration, we show a two-dimensional interface introduced as a series of points which form line segments. To construct the initial boundary, we try to keep the number of points representing the interface to a minimum to avoid long computational time during the calculation of the initial level set function. Figure 5.2 shows an example to illustrate the technique used for the calculation of a common distance function. Consider the grid point $M$ whose distance from the boundary we want to calculate. We first calculate the distance of $M$ from all segments which form the boundary. The line segments are $AB$, $BC$, $CD$, $DE$, and $EF$. From elementary mathematics we know that the distance of a point from a line is the length of the normal line drawn from point to the line, e.g., the distance of $M$ from line segment $CD$ is $MG$. However, there is not always an intersection point between the normal line and segment as it is the case for the line segment $EF$, where $MH$ is outside $EF$. In such cases the distance is the minimum of $ME$ and $MF$.

Figure 5.2: Illustration of the calculation of a common distance function of a grid point from a boundary.

After the calculation of these different distances and getting the minimum value among them, we still need to assign a sign to this value depending on the position of the grid point relative to the boundary. In order to do this we choose a reference point $R$ inside the boundary and at the bottom of the simulation domain. If the connection line between a grid point and the reference point intersects the boundary at an odd or even number of times, the grid point is outside or inside the boundary and the common distance function must be multiplied by plus or minus one, respectively. For example, the grid point $N$ is outside the boundary since the line $RN$ intersects the boundary once at $I$ and the grid point $O$ is inside the boundary, because the line $RO$ intersects the boundary twice, at $J$ and $K$.

Since the level set algorithm will later only work in a narrow band and the calculation of the distance function especially in three dimensions is very CPU expensive, it is sufficient to perform the distance calculations near the initial boundary. This can be achieved, e.g., by a recursive algorithm walking along the boundary. First, a point whose distance to the boundary is smaller than the width of the future narrow band is identified by walking along the boundary of the simulation domain in $y$-direction or $z$-direction for two or three dimensions, respectively. When such a point is found, it is used as the starting point for the recursion. In the recursion only points with a distance smaller than the narrow band width are considered for the next step. The maximum distance ever computed during this procedure is used to initialize the grid points above and below the boundary with a positive or negative sign, respectively. This method reduces the computational effort of initialization from $O(N^3)$ and $O(N^2)$ to $O(N^2)$ and $O(N)$ in three and two dimensions, respectively, where $N$ is the number of grid points in each direction.

For three-dimensional initialization, in analogy to two dimensions one has to calculate the distance of the grid points to different triangles which represent the boundary. The distance of a point from a triangle $t$ is given by $d(P,t)=\min_{T\in t} \Vert P-T \Vert _2$ and the distance $d(P,b)$ of a point from the initial boundary $b=\{t_1,\ldots,t_m\}$, consisting of a set of triangles, is given by $d(P,b)=\min_{t\in b} d(P,t)$.

Figure 5.3: The signed distance function is used as the initial level set function which corresponds to the initialization of a rectangular trench. The grid resolution was 80$\times$160.

The initial level set function for a rectangular trench is shown in Figure 5.3. As it can be seen in this figure the signed distance function is a linear function with a value of zero at the boundary.

5.3 Transport of Particles

Before we begin to describe the radiosity model used by transport of particles, we have to present the definition of the visibility of a point from another point. This plays an important role for the calculation of the flux at the surface elements as it can be seen later in this section.

Figure 5.4 is an illustrative example for the definition of visibility. There are different kinds of visibility which have to be defined. The first kind is the visibility among the points of a boundary such as $(A,B)$, $(A,C)$, and $(B,C)$. The second kind is the visibility between a boundary point and a source point which is assumed to be $S$ in Figure 5.4. Whereas the pairs of points $(A,S)$, $(A,B)$, $(B,S)$, and $(B,C)$ are visible from each other, there is no visibility between the pairs of points $(A,C)$ and $(S,C)$.

Figure 5.4: An illustrative example for the definition of the visibility between two points.

A formulation of the radiosity model [80] for the case of luminescent reflection, pertaining to low energy particles (cf. Section 4.1), can be used in the deposition simulations.

To model deposition it is assumed that the distribution of the particles coming from the source obeys a cosine function around the normal vector of the plane in which the source lies [19]. This implies that the incoming flux at a surface element is proportional to the cosine of the angle between the connecting line between the center of mass of a surface element and the source, and the normal vector of the source plane.

The flux reaching the surface elements obtained by the surface extraction step may be written as a vector,

$$\mathrm{Flux} = \beta_0 I_S + \beta \Psi L I_R$$

$$= \frac{\beta-\beta_0}{1-\beta}I_S +\frac{\beta(1-\beta_0)}{1-\beta} \underbrace{L^{-1} (L^{-1} - (1-\beta)\Psi)^{-1}}_{T:=} I_S.$$

Here $I_S$ is the vector of fluxes coming from the sources to the surface elements, $I_R$ is the vector of fluxes that arrive because of reflections, $\beta_0$ the sticking coefficient for particles coming directly from the source, $\beta$ the sticking coefficient for secondary bounces, $L$ the diagonal matrix containing the lengths of the surface elements, and

$$\Psi_{ij} = \frac{n_i\cdot(t_j-t_i) n_j\cdot(t_i-t_j)}{\pi \vert t_j-t_i\vert^3}v_{ij} %[\mbox{$i$ visible $j$}]$$

where $t_i$ are the centroids of the surface elements, $n_i$ their unit normal vectors, and $v_{ij}$ is $1$ if the surface element $j$ is visible from $i$ or 0 if not.

The second line in the equation above is obtained from the first one and the relationship $I_R= (1-\beta_0)I_S + (1-\beta) \Psi L I_R$ results from straightforward algebraic manipulations. In the case of multiple, low energy species, the calculation of the visibility matrix and the inverse $T$ only depends on topographic information and thus does not have to be repeated for each species.

Most of the computation time for the simulation of the transport of the particles above the wafer by the radiosity model is consumed by first determining the visibility between a surface element and the source, and secondly, between two different surface elements. The latter has an $O(n^{2})$ operation complexity, where $n$ is the number of surface elements growing approximately with $O(N^{2})$. If the connecting line between the center of mass of two surface elements does not intersect the surface, i.e., the zero level set, those surface elements are visible from each other. In order to decrease the computational effort related to determining the visibility between the surface triangles, we have assumed that two triangles are visible from each other if the center point of the grid cells in which the triangles are located, are visible from each other. Since there are at least two triangles in each grid cell, considerable time is saved.

5.4 Extending the Speed Function

In most of applications the speed function is only given on the boundary. For example, in the simulation of etching and deposition processes, the speed function on the boundary depends on the visibility of the source of species from the boundary segments. However, this visibility can only be calculated between the boundary and source of species, since it does not have any physical meaning, when we want to calculate it between the other level sets and the source of species [2,3,4]. Therefore, we must extend the speed function.

After we have calculated the speed function on the boundary, where the models of calculation differ depending on the process conditions (in Chapter 4 we will introduce several different models which are used to calculate the speed function on the boundary), we can start with extending the speed function. An extended speed function must satisfy two basic requirements. First it has to match the speed function calculated from a physical model on the boundary, and second it has to guarantee that it moves the other level sets in such a way that the signed distance function is preserved.

Assume that $\varphi (\mathbf{x},0)$ is a signed distance function which satisfies $\vert\nabla\varphi (\mathbf{x},0)\vert=1$. We will prove that $\vert\nabla\varphi (\mathbf{x},t)\vert$ will be equal to one, i.e., the level set function remains a signed distance function, if we build the following equation to extend the speed function [104]

$$\nabla\varphi \cdot\nabla F_{ext}=0.$$ (5.1)

In order to prove this we build the following equation

$$\frac{d\vert\nabla\varphi \vert^{2}}{dt}=\frac{d(\nabla\varphi \cdot \nabla\varphi )}{dt}=2\nabla\varphi \cdot \frac{d\nabla\varphi }{dt}.$$ (5.2)

The derivation of the level set equation in space is

$$\nabla(\varphi _{t})+\nabla F_{ext}\cdot \vert\nabla \varphi \vert+\nabla\vert\nabla \varphi \vert\cdot F_{ext}=0$$

which is represented as follows:

$$\frac{d\nabla\varphi }{dt}=-\nabla F_{ext}\cdot \vert\nabla\varphi \vert-\nabla\vert\nabla\varphi \vert\cdot F_{ext}.$$ (5.3)

Substituting $d\nabla\varphi /dt$ from (5.3) in (5.2) gives

$$\frac{d\vert\nabla\varphi \vert^{2}}{dt}=-2\nabla\varphi \cdot(F_... ...t \nabla\vert\nabla\varphi \vert+\vert\nabla\varphi \vert\cdot\nabla F_{ext})$$

which can be rewritten to the following form

$$\frac{d\vert\nabla\varphi \vert^{2}}{dt}=-2\nabla\varphi \cdot\na... ...a\varphi \vert-2\nabla\varphi \cdot\nabla\vert\nabla\varphi \vert\cdot F_{ext}.$$ (5.4)

The first term at the right hand side of (5.4) is zero because of (5.1) and if $\vert\nabla\varphi (\mathbf{x},t)\vert$=1, the second term on the right hand side vanishes simultaneously with the term on the left. Therefore, we have proven that the level set function remains a signed distance function for all time if we extract the speed function using (5.1).

Now the speed function can be extended. For a level set function $\varphi _{n}$, where $n$ stands for the $n$th time step, we choose a signed distance function $\varphi _{temp}$ whose zero level set is equal to the zero level set of $\varphi _{n}$. As mentioned above, we need to satisfy the following equation

$$\nabla\varphi _{temp}\cdot \nabla F_{ext}=0.$$

Note that the temporary signed distance function $\varphi _{temp}$ is only used to extend the speed function and to guarantee that the signed distance function is preserved. It plays no other role such as for re-initialization, etc. The calculation of $\varphi _{temp}$ is done by solving the so called Eikonal equation

$$\vert\nabla\varphi _{temp}\vert=1$$

which is solved very effectively using the fast marching method. In the next section we will describe this method shortly.

5.5 Fast Marching Method

The fast marching method [78] is a method for the efficient calculation of the signed distance function. As mentioned in Section 2.1.3, the fast marching method is used in the boundary value formulation where the speed function must have a unique sign. Therefore, it is necessary to calculate once the signed distance function above the boundary and once below the boundary. The goal is to solve the Eikonal equation $\vert\nabla T\vert F=1$ using proper and efficient algorithms. The key idea is to build an approximation of the gradient term which correctly deals with the development of corners and cusps in the evolving solution. It is well-known that the Eikonal equation becomes non-differentiable, and an appropriate weak solution must be built which is related to the entropy condition of propagating interfaces introduced in [76]. One of the simplest entropy-satisfying approximation of the gradient is from Godunov, and was used, for example, by Rouy and Tourin [68] to solve the Eikonal equation as follows:

$$\sqrt{\max(D^{-x}_{ijk}T,-D^{+x}_{ijk}T,0)^{2}+ \max(D^{-y}_{ijk}... ...{+y}_{ijk}T,0)^{2}+ \max(D^{-z}_{ijk}T,-D^{+z}_{ijk}T,0)^{2}}=\frac{1}{F_{ijk}}$$ (5.5)

where

$$D^{-x}_{ijk}T=\frac{T(i,j,k)-T(i-1,j,k)}{\Delta x}$$

$$D^{+x}_{ijk}T=\frac{T(i+1,j,k)-T(i,j,k)}{\Delta x}$$

and $i$, $j$, and $k$ are the indices of a grid point on the $x$, $y$, and $z$ axis, respectively. $D^{-y}$, $D^{+y}$, $D^{-z}$, and $D^{+z}$ are defined in analogy to $D^{-x}$ and $D^{+x}$. Additional schemes for solving Hamilton-Jacobi equations may be found in [11,63].

The upwind scheme^5.1 of (5.5) enables the calculation of $T$ from the grid points adjacent to the boundary which have the smallest values of $T$. The boundary is swept along considering points in a narrow band (cf. Section 5.6) around the boundary. Marching this narrow band forward and freezing the values of the existing points and bringing the new points into the narrow band results in a solution of the Eikonal equation within the narrow band. More details can be found in [78,79].

5.6 Narrow Banding Technique

The straightforward method for solving the level set equation is to solve it in the entire computational domain where one has to update all the level sets and not only the zero level set. This approach has been called the full matrix method. The advantage of this simple method is the simplicity of data structures and operations. Furthermore, it is a good starting point for the construction of a level set code. There are even cases in which this type of implementation is necessary, such as if all the level sets are themselves important as is the case in problems encountered in image processing. However, there is the disadvantage that an increase in the resolution of the grid increases the calculation time and the complexity according to $O(N^{2})$ and $O(N^{3})$, in two and three dimensions, respectively.

One technique to reduce the calculation time is to use adaptive mesh refinement. This refinement may be required especially in two regions. The first region is where the curvature of the boundary is high and the second is where the speed function changes very rapidly. For example, if the zero level set identified with the boundary is the object of interest, as is normally the case, then one has to adaptively refine the mesh around the location. However, considerable care must be taken at the interfaces between the fine and coarse cells. In particular, a subtle update strategy for the level set function values is required at so called hanging nodes where the boundary between two levels of refinement does not have a full set of nearest neighbors. Whereas advection terms which do not depend on the curvature lead to hyperbolic equations and a straightforward interpolation of the level set values from the coarse grid points easily produces the level set values at hanging nodes, there are some complications for the curvature dependent advection terms. The situation is not so straightforward, since this corresponds to a parabolic term that can not be approximated by simple interpolation [80].

The idea leading to fast level set algorithms stems from the observation that only the values of the level set function near its zero level set are essential, and thus only the values at the grid points in a narrow band around the zero level set have to be calculated. The advantages of this technique are as follows:

• Efficiency: when working only in a narrow band around the zero level set, the operation complexity reduces from $O(N^{3})$ and $O(N^{2})$ to $O(lN^{2})$ and $O(lN^{})$, for three and two dimensions, respectively. Here, $l$ is the number of grid points in the narrow band.
• Extending the speed function: as mentioned in Section 5.4, in applications linking to physical models the speed function is not known on the whole simulation domain, but only at the surface. In order to use the level set method it has to be suitably extended from the known values to the whole simulation domain. Therefore, the problem of obtaining a smooth speed function within the entire computational domain is avoided or reduced into the narrow band.
• Stability: one of the most important things which must be taken into account, is the stability problem during the update of the level set function using finite difference methods where the CFL (cf. Section 5.13) condition must be fulfilled. Whereas this condition can be satisfied more easily by the zero level set and adjacent level sets, it is difficult to fulfill it for all level sets.

Figure 5.5: Simulation of void formation for a deposition process. A level set grid of $80\times 160$ points was used.

Figure 5.5 shows the simulation results of different steps of a deposition process into a typical trench. Except for step 0 for which the level set function with and without narrow banding has been shown, other intermediate level set functions corresponding to Figure 5.5 have been shown only within the narrow band in Figure 5.6. The grid resolution is $80\times 160$. The zero level set can be seen easily in each intermediate step and the formation of a void in step 48, as well.

Figure 5.6: The level set functions at step 0, $12$, $24$, $36$, and $48$ during the simulation as shown in Figure 5.5. Inside the narrow band the level set function is retained to the end of the simulation, whereas the level set function values of the other points have been substituted with the width of the narrow band multiplied by $1$ or $-1$ depending on their position.

[Step 0 without narrow banding] [Step 0] [Step $12$] [Step $24$] [Step $36$] [Step $48$]

5.7 Combining Narrow Banding and Extension of the Speed Function

In this section we present the implementation of an algorithm which combines narrow banding and extension of the speed function. This algorithm works as follows: first the initial points near the zero level set, where the speed function is known, and the neighboring trial points are determined. It is checked in the main loop, if there still is a trial point to be considered in the narrow band. All trial points are stored on a heap ordered by their distance from the zero level set. If there is a point to be considered, its distance is approximated, its extended speed is calculated, and its neighbors are updated accordingly. Finally, after the main loop, the bookkeeping information for the narrow band points is updated using the distance information just computed.

The detailed steps using several auxiliary functions are as follows: concerning the data structures, the information about the level set grid, the distance function, the extended speed function, and the tags for the fast marching method are stored in arrays. The trial points are stored on a heap. The implementation flow is as follows:

1. First find the grid points whose speed function values are initially known. These values are computed in the physical simulation step and translated to the grid. Next compute the distance for the initial points, tag them as known, and initialize the corresponding grid points of the speed function. Find the trial points which are the neighbors of the initially known points, and compute their tentative distance values. All other points are far points.
2. While there are trial points, do the following:
   1. Remove the first point from the heap and call it $a$. It has the smallest distance from the zero level set of all points on the heap.
   2. If narrow banding is used and the maximum width of the narrow band is smaller than the distance of $a$, return from the loop.
   3. Mark $a$ as known.
   4. For all neighbors $b$ of $a$, do the following:
      1. If $b$ is a far point, recompute its distance and speed function values and mark it as a trial point.
      2. If $b$ is a trial point, recompute its tentative distance and speed function values, unless it was computed in the previous step.
3. If narrow banding is used, set the distance values of all points which are not marked as known, to the width of the narrow band. Set the sign of the distance function.
4. Finally, return two objects, namely the new signed distance function and the extended speed function.

5.8 Coarsening Algorithm

When using the radiosity model to simulate the transport of particles above the wafer, two operations consume most of the computation time. The first operation is the determination of the visibility between all surface elements, which requires $\binom{n}{2}$ visibility tests, where $n$ denotes the number of surface elements extracted from the level set grid. The second operation is the solution of a certain system of linear equations, which leads to the calculation of the inverse of a dense $n^2\times n^2$ matrix.

Hence it is mandatory to keep the number of surface elements as low as possible. However, decreasing the number of surface elements must be done in such a way that the high resolution is obtained where it is needed, e.g., near the trench opening and at the bottom of the trench. One approach is to devise a refinement and coarsening strategy for unstructured grids at the level set implementation and the algorithms working on it. This, however, complicates the implementation because of the complex fast marching algorithm necessary at the unstructured grid compared to the algorithm at the rectangular grid which is used normally for the implementation of the level set method. In this work a different approach was taken by coarsening the surfaces after they have been extracted from the level set grid.

The coarsening algorithm works by walking down the list of surface elements extracted as the zero level set and calculating the angle between two neighboring surface elements. Whenever this angle is below a certain threshold value of a few degrees, the neighboring elements are coalesced into one. After one sweep through the list, the algorithm can be reapplied for further coarsening. After $k$ coarsening sweeps, at most $2^k$ surface elements are coalesced into one. The resulting longer surface elements are used for the radiosity calculation, after which the fluxes are translated back from the coarsened elements to the original ones. This is a heuristic approach and has shown excellent results.

5.9 Advancing the Level Set Function Using Finite Difference Schemes

To discretize the level set equation, one can substitute the time and spatial derivatives with forward and central differences, respectively. We now consider a boundary with a right angle at the corner that is moved with a constant speed function normal to the boundary. Studies in [80,35] have shown that using the central differences for spatial derivatives leads to false values of the gradient at the corner point. Since the slope $\nabla_{\mathbf{x}}\varphi$ is not defined at corners, the central difference approximation sets it to the average of the left and right slopes. Therefore, this wrong calculation of the slope propagates away from the corner and leads to oscillations.

An alternative approach would be to add a viscosity term to the right hand side of the level set equation as follows [80]:

$$\varphi _{t}+F(t,\mathbf{x})\vert\vert\nabla_{\mathbf{x}}\varphi \vert\vert=\varepsilon \Delta\varphi$$

where $\varepsilon$ is a positive constant. Although the solutions $\varphi _{\varepsilon }$ guarantee the smoothness and their limits yield an appropriate weak solution as $\varepsilon \to 0$, the smoothness introduced by $\varepsilon \Delta\varphi$ often results in significant rounding of corners.

The best approach to discretize the equations of the form $\varphi _{t}+(G(\varphi ))_{\mathbf{x}}$ are the following space convex schemes, where it is assumed that the flux $G(\varphi )$ is convex, i.e., $d^{2}G/d\varphi ^{2}>0$ [80]. These approaches ensure that discontinuities and boundaries remains sharp.

The simplest scheme to discretize the level set equation is the first order space convex scheme which is as follows:

$$\varphi ^{n+1}(t)=\varphi ^{n}(t)-\Delta t[\max(F_{ijk},0)\nabla^{+}+\min(F_{ijk},0)\nabla^{-}]$$ (5.6)

where the $\varphi ^{n+1}(t)$ and $\varphi ^{n}(t)$ are the level set functions at the $(n+1)$th and $n$th time steps, respectively, and

$$\nabla^{+}= \Bigl[\max(D^{-x}_{ijk}\varphi ^{n},0)^{2}+\min(D^{+x}_{ijk}\varp... ...^{2}+\max(D^{-y}_{ijk}\varphi ^{n},0)^{2}+\min(D^{+y}_{ijk}\varphi ^{n},0)^{2}+$$

$$\max(D^{-z}_{ijk}\varphi ^{n},0)^{2}+\min(D^{+z}_{ijk}\varphi ^{n},0)^{2}\Bigl]^{\frac{1}{2}}$$

$$\nabla^{-}= \Bigl[\max(D^{+x}_{ijk}\varphi ^{n},0)^{2}+\min(D^{-x}_{ijk}\varp... ...^{2}+\max(D^{+y}_{ijk}\varphi ^{n},0)^{2}+\min(D^{-y}_{ijk}\varphi ^{n},0)^{2}+$$

$$\max(D^{+z}_{ijk}\varphi ^{n},0)^{2}+ \min(D^{-z}_{ijk}\varphi ^{n},0)^{2}\Bigl]^{\frac{1}{2}}.$$

For the second order space convex scheme $\nabla^{+}$ and $\nabla^{-}$ are defined as follows:

$$\nabla^{+}=\sqrt{\max( A ,0)^2 + \min( B ,0)^2 +\max( C ,0)^2 + \min( D ,0)^2 +\max( E ,0)^2 + \min( F ,0)^2}$$

$$\nabla^{-}=\sqrt{\max( B ,0)^2 + \min( A ,0)^2 +\max( D ,0)^2 + \min( C ,0)^2 +\max( F ,0)^2 + \min( E ,0)^2}$$

where

$$A := D_{ijk}^{-x}\varphi ^{n} + \frac{\Delta x}{2} m( D_{ijk}^{-x-x}\varphi ^{n}, D_{ijk}^{+x-x}\varphi ^{n})$$

$$B := D_{ijk}^{+x}\varphi ^{n} - \frac{\Delta x}{2} m( D_{ijk}^{+x+x}\varphi ^{n}, D_{ijk}^{+x-x}\varphi ^{n})$$

$$C := D_{ijk}^{-y}\varphi ^{n} + \frac{\Delta y}{2} m( D_{ijk}^{-y-y}\varphi ^{n}, D_{ijk}^{+y-y}\varphi ^{n})$$

$$D := D_{ijk}^{+y}\varphi ^{n} - \frac{\Delta y}{2} m( D_{ijk}^{+y+y}\varphi ^{n}, D_{ijk}^{+y-y}\varphi ^{n})$$

$$E := D_{ijk}^{-z}\varphi ^{n} + \frac{\Delta z}{2} m( D_{ijk}^{-z-z}\varphi ^{n}, D_{ijk}^{+z-z}\varphi ^{n})$$

$$F := D_{ijk}^{+z}\varphi ^{n} - \frac{\Delta z}{2} m( D_{ijk}^{+z+z}\varphi ^{n}, D_{ijk}^{+z-z}\varphi ^{n}).$$

The switch function $m$ is defined as

$$m(a,b):= \begin{cases} \begin{cases} x & \text{if $\vert x\vert\l... ...t$} \end{cases}& \text{if $xy\ge0$},\\ \quad 0 & \text{if $xy<0$}. \end{cases}$$

This scheme takes care of shocks by building an appropriate switch. A comparison between these two different schemes can be found in [35].

5.10 Iterative Solver

As mentioned in Section 5.3 the radiosity model assumes that the total flux depends on the flux directly from the source, as well as an additional flux due to the particles which do not stick and are re-emitted. After discretizing the problem the flux vector whose elements are the total flux at different surface elements can be expressed by a matrix equation.

There are two numerical approaches to solve this equation. The first is to use a direct solver for the matrix equation. Although this is practical in two dimensions [13,14,11], it becomes impractical due to the computational effort needed to calculate the inverse matrix for three-dimensional problems. In three dimensions the equation is solved iteratively.

The iterative solution developed by Adalsteinsson and Sethian [80] consists of a series expansion of the radiosity matrix. Suitably interpreted, it can be viewed as a multi-bounce model in which the number of terms in the series expansion corresponds to the number of bounces that a particle can undergo before its effect is negligible.

The iterative solution allows one to check the error term to determine how many terms must be kept. Since most of the particles either stick or leave the domain after a reasonable number of bounces, this is an efficient approach.

The reflected intensity after the $k$-th bounce is defined as $I_{R,k}$. The relationship between reflected intensity and source intensity at the $i$-th surface element is written (cf. Section 5.3) as follows:

$$I_{R,1}^{i}=(1-\beta_0)I_{S}^{i}.$$ (5.7)

In matrix form, this becomes

$$I_{R,1}=(1-\beta_0)I_{S}$$ (5.8)

$$I_{R,k+1}=(1-\beta)\Omega I_{R,k}$$ (5.9)

where $\Omega$ is a matrix whose elements stand for the influence of the visibility term, diffuse reflection, effective surface area of a surface element, and the distance between the center of mass of two surface elements [80].

Now $I_{S,k}$ is defined to be the portion that sticks after the $k$th bounce.

$$I_{S,0}=\beta_0 I_{S}$$ (5.10)

$$I_{S,1}=\beta \Omega( 1-\beta_0)I_{S}.$$ (5.11)

In general,

$$I_{S,k+1}=\frac{\beta}{1-\beta}I_{R,k+1}=(1-\beta)\Omega I_{S,k}.$$ (5.12)

Therefore, by going back to the initial expression (5.10)

$$I_{S,k}=\beta(1-\beta)^{k-1}(1-\beta_0)\Omega^{k}I_{S}$$ (5.13)

and thus the total intensity after $N$ reflections is given by

$$I_{N}=\beta(1-\beta_0)\Bigl[\sum_{k=1}^{N}(1-\beta)^{k-1}\Omega^{k}\Bigl]I_{S} + \beta_0 I_{S}.$$ (5.14)

Each application of the operator may be viewed as either an additional term in the expansion or an additional included bounce. It is important to note that there is a recurrence relation for $I_{N}$ given by

$$I_{N+1}=\beta(1-\beta_0)\Bigl[\sum_{k=1}^{N+1}(1-\beta)^{k-1}\Omega^{k}\Bigl]I_{S}+\beta_0 I_{S}$$ (5.15)

and then

$$I_{N+1}=\beta(1-\beta_0)\Bigl\{(1-\beta)\Omega\Bigl[\sum_{k=1}^{N}(1-\beta)^{k-1}\Omega^{k}\Bigl]+\Omega\Bigl\}I_{S}+\beta_0 I_{S}$$ (5.16)

$$I_{N+1}=(1-\beta)\Omega(I_{N}-\beta_0 I_{S})+\beta(1-\beta_0)\Omega I_{S}+\beta_0 I_{S}$$ (5.17)

$$I_{N+1}=(1-\beta)\Omega I_{N} +(\beta-\beta_0)\Omega I_{S} +\beta_0 I_{S}.$$ (5.18)

By constructing the remainder term $I_{N+1}-I_{N}$, the convergence of the expansion can be measured and sufficient terms can be kept to bound the error below a user-specified tolerance.

5.11 Input File of Three-Dimensional ELSA

Three-dimensional ELSA requires a library called IPD (input deck library) developed by the Institute for Microelectronics [51]. It is used for parsing input files and must be installed first. An example of an IPD input file including parameters of the simulator, process parameters, and the initial geometry is as follows:

Simulation
{
  Name = "Sample";
  Comment = "a sample deposition with void";
  OutputDirectory = "./output/";
  WriteAfterEachIteration = "yes";
  PrintConfig = "short";
}

// Parameters of the level-set algorithm
LevelsetAlgorithm
{
  nx=30;
  ny=30;
  nz=30;
}

Sources
{
  SourceGroup1
  {
    z = 2.0;
    x1 = 0.0;
    y1 = 0.0;
    x2 = 1.0;
    y2 = 1.0;
    xSourceCount = 4;
    ySourceCount = 4;
    Intensity = 20;
    n = 1;
    Direction = [0,0,-1];
  }
}

Steps
{
  Deposition1
 {
    Name = "Deposition1";
    Material = "Si3N4";
    Time = 8.0;
    MaxThickness = 0.23;
    SourceStickingProbability = 0.1;
    ReflectionStickingProbability = 0.1;
  }
}

// We define the initial geometry here as a list of triangles:
InitialBoundary
{
  lowerLeft= [0.0, 0.0, 0.0];
  upperRight=[1.0, 1.0, 1.0];
  Points
  {
    wafer\_level = 0.75;
    p11=[~InitialBoundary.lowerLeft[0],
         ~InitialBoundary.lowerLeft[1],
         wafer\_level];
    p12=[~InitialBoundary.lowerLeft[0],
         ~InitialBoundary.upperRight[1],
         wafer\_level];
    p13=[~InitialBoundary.upperRight[0],
         ~InitialBoundary.upperRight[1],
         wafer\_level];
    p14=[~InitialBoundary.upperRight[0],
         ~InitialBoundary.lowerLeft[1],
         wafer\_level];
    p21=[0.45, 0.45, wafer\_level];
    p22=[0.45, 0.55, wafer\_level];
    p23=[0.55, 0.55, wafer\_level];
    p24=[0.55, 0.45, wafer\_level];
    p31=[0.4, 0.4, 0.2];
    p32=[0.4, 0.6, 0.2];
    p33=[0.6, 0.6, 0.2];
    p34=[0.6, 0.4, 0.2];
   }
  Boundary
   {
    q1=["p11","p12","p22","p21"]; //q1 ... q4 top planes
    q2=["p12","p13","p23","p22"];
    q3=["p13","p14","p24","p23"];
    q4=["p14","p11","p21","p24"];
    q5=["p21","p22","p32","p31"]; //q5 ... q8 side planes
    q6=["p22","p23","p33","p32"];
    q7=["p23","p24","p34","p33"];
    q8=["p24","p21","p31","p34"];
    q9=["p31","p32","p33","p34"]; //q9 bottom plane
  }
}

In the next section we present another program library suitable for data exchange between different process simulators even when they are based on different native file formats.

5.12 Fully Three-Dimensional Process Simulation

ELSA is a stand alone application for three-dimensional topography simulation and uses IPD for input file as mentioned in previous section. Furthermore, TOPO3D has been also developed. The kernel of TOPO3D is based on ELSA, but is also linked to a program library, that handles objects for fully three-dimensional semiconductor process simulations. This program library is called WSS (wafer state server) [14].

The idea of WSS stems from necessity that the simulation informations generated by a simulator must be usable not only by the same simulator but also by other tools and simulators. These tools can be visualization programs, for instance. Mostly the concrete syntax based on the informations are saved, is not important for users. However, a low memory consumption and a fast access during writing and reading an input or output file are very essential. Different simulators save the simulation informations in different file formats which are not compatible to each other. Therefore, in general a certain file format can only be read by the related tool. This incompatibility can be overcome with introducing the common data format presented by WSS [45].

WSS is a program library and file format for handling three-dimensional objects in semiconductor process simulations developed at our institute [14]. It is a solution for the integrated simulation of three-dimensional manufacturing processes. A generic data model suitable for process and device simulations allows an efficient data exchange between simulators even when they are based on different native file formats. It is also able to handle different meshes and distributed quantities stored thereon. The program library also defines algorithms to perform geometrical operations for fully three-dimensional process simulations as they are used in topography simulations.

Figure 5.7: Overview of the simulation flow in TOPO3D.

Figure 5.7 shows an overview of the simulation flow achieved when using TOPO3D. The simulation begins at the WSS. WSS provides us with the initial surface which is then propagated. The initial surface is represented by a set of triangles. In the next step the level set function is initialized with the signed distance function. After this step the time stepping is started. At each time step a data exchange between the simulation model and the level set kernel takes place. The simulation model provides the speed function needed by the level set kernel to propagate the surface. Vice versa the model also needs information about the actual distance function from the surface propagation step. After $n$ time steps the final surface is extracted as the zero level set of the level set function. After coarsening the final surface and merging it with the initial surface, a meshing step is performed and the information is stored in WSS.

5.13 Stability of the Simulator

To advance the level set function we have used a second order space convex finite difference scheme [24,20,19,13]. Consider $\Delta x$, $\Delta y$, $\Delta z$, and $\Delta t$ as discretization steps in space and time, respectively. As mentioned in Section 5.6, a necessary condition for the stability of this scheme is the CFL condition which requires that

$$\Delta t \cdot F_{\max}\leq \min(\Delta x,\Delta y, \Delta z)$$

where $F_{\max}$ is the maximum value of $F$ at the grid points. The CFL condition guarantees that the front can not cross more than one grid cell during each time step. In order to have a stable simulator based on the finite difference method, the CFL condition must be satisfied [69].

However, there is a problem stemming from the CFL condition, which limits the simulator performance. If we increase the spatial resolution by $\lambda$, assuming that $F_{\max}$ remains constant, we have to reduce the maximum $\Delta t$ by the same factor $\lambda$, which increases the number of simulation steps by $\lambda$ to reach the same thickness. Furthermore, an increase in spatial resolution by $\lambda$ approximately increases the number of extracted surface elements by $\lambda^{2}$ in three dimensions. Therefore, the computational effort of the visibility determination is increased by $\lambda^{4}$. In total, an increase in spatial resolution by $\lambda$ leads to a minimal increase of simulation time by a factor $\lambda^{5}$, if the most precise visibility calculation is used.

5.14 Summary

State of the art algorithms for surface evolution processes like deposition and etching processes in two and three dimensions have been implemented. A general purpose topography simulator has been developed based on the level set method combining the narrow banding and fast marching methods to extend the speed function. The speed of the simulation has been improved in several steps, e.g., in initialization, visibility determination, and solving the radiosity matrix. To be able to handle the objects for fully three-dimensional semiconductor processes, TOPO3D which has ELSA as its kernel has been developed and is linked to the WSS program library.

Footnotes

... scheme^5.1

An upwind scheme correctly respects the upwind nature of the differential equation and sends information in the direction that correctly matches the differential equation. Furthermore, an associated issue is the stability of a scheme. Stability is examined by finding those values for the time step and space step such that the small errors in the solution do not grow uncontrollably. Typically, the stability of a scheme depends on a balance between the time step $\Delta t$, the space step $\Delta \mathbf{x}$, and the speed; this is known as the CFL (Courant-Friedrichs-Levy) condition. For a constant speed equal to one, stability of the upwind scheme will require that $\Delta t\le \Delta \mathbf{x}$.