Step flow-based cellular automaton for the simulation of anisotropic etching of complex MEMS structures

An enhanced octree data representation is developed for use in atomistic simulations of the evolution of complex multivalued surfaces appearing during anisotropic etching of crystalline silicon in micro electro mechanical systems (MEMS) applications. The octree provides a good balance between accuracy, memory efficiency and calculation speed. In combination with the octree, a step flow-based cellular automaton (CA) model is considered, which can be used to convert the experimental macroscopic etch rates into atomistic removal rates for direct use in the simulations. This involves minimizing the global error between the experimental values and the non-linear analytical expressions of the etch rates for a small set of chosen orientations. The orientations where the etch rate suffers a sudden change against the cutting angle play a crucial role to generate an optimal solution. The simulated etch rates for KOH and KOH/IPA systems in different concentrations show good agreement with the experiments. Due to the improved description of the anisotropy, the propagation of the etch front in realistic engineering applications can be simulated accurately.


Introduction
Etching represents a flexible micromachining technique, with many different choices of etchant and etching conditions to obtain a desired structure. Combined with other micromachining methods, it enables the fabrication of a wide variety of micro electro mechanical systems (MEMS) with diverse functionalities. Sometimes the goal is to avoid the undercutting effects caused by the fast, high Miller index orientations which result from the etch anisotropy. However, sometimes those same planes can be utilized to create many complex threedimensional (3D) structures with intricate details [1,2].
At present, many tools and models for the simulation of etching have been developed [3]- [5]. Although these methods are powerful, they still present many difficulties in practice. In the case of the atomistic methods, such as the cellular automata (CAs) and the kinetic Monte Carlo (KMC) simulators, fulfilling the topological requirements for the data representation of a 3D structure during the dynamic evolution of the surface is still an important issue. The surface advancement and the cell removal model used in current simulators [6,7] are time-consuming and require large use of memory. Without an efficient model, some essential properties of the surface structure cannot be described in the simulations. In [8], the propagation of the etch front is described by a Hamiliton-Jacobi equation combined with an etch rate interpolation model based on values of the (100), (110) and (111) orientations. In this study, the presented step flow CA model provides more flexibility to simulate various enchants in different concentrations. In addition, the atomistic solution is computationally less intense than the solution of the Hamiliton-Jacobi equation.
Another issue is the conversion of the macroscopic anisotropy of the etch rate into the calculation principles and parameters considered in an atomistic model. The current CA models provide a rather limited analysis of the atomistic structure of the etched surfaces [3]- [5]. In the 3 best case, some models consider a few orientations such as (100), (110), (111) and (311) or (411) in order to operate all types of atoms. Although some models have the ability to fit 5 or 6 orientations [5], they need a complex trial and error process to calibrate the final parameters. This type of limitation leads, in general, to both poor fits for the development of new etchants and bad simulation accuracy.
Additionally, the fabrication of micro devices is based on composite micromachining processes which typically combine the use of anisotropic etching, multi-step mask transfers, thin film growth and deposition, and/or deep reactive ion etching (DRIE). However, most of the available simulation tools can only run these processes through the use of independent modules. A versatile simulator needs a good underlying model that supports the simulation of complex, integrated micro machining processes.
Based on geometrical and kinetic aspects of anisotropic etching of silicon as a step flow process on any surface orientation, we have recently presented the analytical solution of the continuous cellular automaton (CCA) [9]. This is a generalized CA where the occupation of the substrate sites is regarded as a continuous variable in the range [0, 1], in contrast to discrete occupation, where only two states are possible, 1 (occupied) and 0 (empty) [3]. Gosálvez et al [9] focuses mainly on: (i) the presentation of a new classification of the surface sites and the atoms occupying them, enabling the treatment of step flow as a mean field process involving alternating removals of complete rows of atoms, (ii) the derivation of the analytical expressions for the etch rates of the different surface orientations as a function of the removal rates {r i } of the classified atoms, i.e. R {hkl} = F({r i }), and (iii) the determination of analytical and/or numerical expressions for the removal rates {r i } as functions of the etch rates {R {hkl} }, allowing calibration of the atomistic rates according to the experimentally measured etch rates for a group of surfaces. For simplicity, we refer to this step flow-based CCA as the 'step flow CA model'. The step flow CA model can be viewed as an improved CCA method based on the assumption of step flow propagation of the etching front. One of its most significant features is an improved ease of calibration, enabling a comprehensive transformation of the experimental macroscopic etch rates into atomistic removal rates.
In the present report, we consider the previous data storage and versatility requirements as important aspects for the implementation of any KMC and CA method, and we use the step flow CA model as a specific example to show how these considerations are put in practice. In particular, we focus on the presentation of: (a) an enhanced octree data structure that simultaneously satisfies (a.1) the topological requirements for representing a substrate with an evolving surface using minimal memory, and (a.2) the versatility requirements for simulating composite micromachining processes involving different types of growth, wet and dry etching, as well as the deposition of functional and/or sacrificial layers and masks; and (b) an alternative calibration method for the atomistic rates of the step flow CA model, where we minimize the global error between the experimental and the simulated etch rates. In particular, (b) provides a rather straightforward way to convert the macroscopic anisotropy of the etch rate into values for the parameters used in the simulations. In this manner, the paper reports two aspects of the step flow CA model which are not considered in [9].
Although the paper is not intended as a full presentation of the step flow CA model, which has been reported in detail in [9], the paper does provide a useful general overview of the model, introducing the main concepts and features. We feel that this is particularly beneficial for readers who lack the time and/or patience to browse through the many details of [9], making the step flow CA model accessible to them.

4
The paper is organized as follows. In section 2, we describe different options for data storage leading to the choice of the octree, which is described in detail, including a number of enhancements. In order to appreciate the value of some of the added features, the basic structure of a CA is described in this section. It is however in section 3 that the step flow CA model is reviewed, including the practical aspects of the numerical implementation. This is followed by the presentation of the new method for the calibration of the atomistic removal rates of the step flow CA model in section 4. The results of two multi-process simulations that use the octree and the new calibration method are presented in section 5. Finally, we present the conclusions of the study in section 6.

Data representation
The representation scheme used in a system determines the efficiency of the data storage, the computing costs and the modeling abilities. In order to describe the dynamic movement of the etched surface and the displacement of silicon atoms, two conventional methods have often been used in previous studies. One of them is the so-called dynamic method [3]. Instead of allocating memory for the whole silicon substrate using a 3D array, the method allocates and de-allocates the memory on a layer-by-layer basis, as required by the system evolution. This method is very useful for the case of vertically propagating surfaces, where allocation of new layers always occurs at the bottom region and de-allocation of already used, empty layers occurs at the top region. However, if the simulation uses masking patterns, these effectively pin the top layer of atoms, resulting in an increase of memory with time as the bottom surface propagates vertically and more layers are allocated while none is de-allocated at the top.
Another conventional way for memory management is the use of a list of surface atoms, where new (old) atoms are inserted (removed) during the simulation [5]. This method is useful for atomistic models where the process rates do not depend on the state of the neighborhood or, if such dependence exists, the neighborhood is very small, e.g. one or two atoms. If the neighborhood of an atom is larger, many searches and maintenance work need to be realized inside the list of surface atoms in order to read/write the rates and other parameters required by the model, making the simulations very slow.
As an alternative to the previous methods, there are algorithms which can avoid the representation of empty portions of the system far from the surface (i.e. in the bulk and in the etchant) as well as provide fast simulations even in the case of larger neighborhoods. Two examples are the sparse matrix and the octree structure [10]. For the purpose of simulating complex MEMS engineering applications, the etch front resulting from etching is typically a multi-valued surface, not a solid-on-solid system. This means that the surface advances towards arbitrary directions depending on the mask patterns and the anisotropy of the process, often giving rise to overhanging features, bridges and/or negative slopes. In this context, the octree structure has the ability to model the dynamic movement using minimal computing resources. In this paper, we introduce an enhanced, octree-based CA suitable for crystalline silicon.
As a hierarchical representation, the octree structure is a tree in which every branch contains eight sub-branches (or octants), whose status can be one of 'empty', 'full' or 'partially full'. The tree is generated by recursively subdividing the cubic space into octants, effectively splitting the outermost cubic region into 2 n × 2 n × 2 n unit cubes (or leaves), where n is the resolution parameter (or number of levels, i.e. the number of subdivisions) and 2 n is the 1D size of the tree. The root node of the tree represents the entire space and the last level contains the leaves of the tree. In our implementation, each leaf is topologically a unit cube. In practice, each leaf corresponds to an orthorhombic crystallographic unit cell containing a specific collection of atoms (the crystallographic basis) whose number and positions depend on the crystallographic orientation of the surface. In this paper, each unit cell at the leaves level is also referred to as a voxel. Searching an atom in the tree can be implemented as a tree traversal procedure. The computational cost of the octree search is O(log 8 N ), where N is the number of atoms, as compared to O(N ) for a linear search (i.e. when a linear array holds the data instead of a tree) and O(log 2 N ) for a binary tree search. By displaying the octree as a grid, figure 1(a) shows the structure of an octree in action during a simulation. Note the rather empty areas in between the micro needles as well as in the lower region below the exposed surface, where the use of memory has been reduced substantially. With this approach, only the evolving surface and the very adjacent regions are kept in memory with full details while the other regions are stored in a rather sparse representation. Figure 1(b) shows the detailed contents of the octree leaves (i.e. the atoms) in a small region close to the foot of a micro needle from figure 1(a). For this particular application, we can observe the formation of high-index crystallographic planes such as (311) due to etching.

CCA
CAs are discrete models studied in computational science consisting on a regular grid of cells with certain states. The state of a cell at time t is a function of the states of a finite number of cells (i.e. its neighborhood) at time t − 1. Every cell has the same rule for updating, based on the values in this neighborhood. In this study, we use the so-called CCA, which is a generalization introduced by Zhu and Liu [3] for the CA method. In the CCA method, the silicon substrate is described as a collection of sites that can be fully or partially occupied by the atoms and the propagation of the surface is described as the overall result of gradually decreasing the occupation of the sites according to the site-dependent removal rates of the atoms.
Due to the crystallographic structure of silicon, corresponding to that of diamond (facecentered-cubic lattice with two atoms in the crystallographic basis, one at (0,0,0) and the other 6 at (1/4, 1/4, 1/4)), the minimum element in our substrate is a site (or an atom located in it), not a cell. We restrict the use of the word 'cell' to denote a crystallographic unit cell. We argue that the term 'site' is more general than the term 'atom' since a 'site' may be occupied or empty while the use of 'atom' assumes the occupied state of the site.
The CCA can be defined as a 4-tuple, A = (C, W, F, ), where: 1. C is the CA space. This is the complete collection of all the silicon sites, which may be empty or occupied: Here, n x , n y and n z are the XYZ coordinates of the lattice node that holds the crystallographic unit cell where the site is located, n is the number of levels of the tree, n l is the index of the site inside the crystallographic unit cell and m is the number of sites in the cell. For a substrate whose orientation is given by high Miller indices, m can have a large value (e.g. hundreds). 2. W is the table of nearest neighbors. Due to the crystallographic structure of the silicon crystal, each atom has four neighbors. Thus, W is composed of m matrices of size 4 × 4 in octree space, storing the relative values of the four coordinates of the four neighbors with respect to each atom in the unit cell. For a substrate with a different orientation, W is different and is calculated before constructing the tree. For the n l atom in the cell, the 4 × 4 block in W is: , so that the four neighbors of the target atom can be written as: 3. F is a function that determines the removal rate of the atoms at time t depending on the configuration of their neighborhoods at time t − 1. This can be a lookup table containing the different configurations and their rates, and/or an explicit function of the state of the neighborhood. 4.
is the occupation state of all the sites, which can take any real value in the interval [0, 1], where 0 means 'empty' and 1 means 'full'. Since any real value is possible, partial occupation of the sites is possible. The occupation state of site n x , n y , n z , n l at time t is defined as: π t n x ,n y ,n z ,n l = π(t) = 1 − t 0 r (t)dt, where r (t) is the removal rate of the atom at time t. The removal rate of each atom changes with time according to the changes in the state of the neighborhood. The initial state of π t n x ,n y ,n z ,n l is one, which represents the full occupation of a surface site with an atom newly exposed to the etchant.

Enhancements to the octree
In comparison to a generic tree-based algorithm, our implementation provides some enhancements in order to allow the use of the octree for the simulation of an evolving etched surface. In particular, the tree has special functionalities in order to operate the atoms stored inside it.

Locating the atoms.
Putting in practice the atomistic model based on the step flow aspects of etching (section 3) requires the use of a surface atom classification scheme based on counting the free bonds of the surface atoms. This operation requires a fast locating algorithm for positioning the target atom and its neighbors within the tree. Since this search for neighbors happens for every atom removal, efficient addressability is the key factor to run the hierarchical structure lightly. In our implementation, the algorithm maintains in memory the neighbors of each surface atom (the closest set) so that updating simply involves finding the neighbors under the coherence assumption (i.e. under the assumption that the neighbors are stored in the structure and thus should always be found). The expected time complexity to update the closest set is a constant. The main advantage of this scheme is that the site coordinates are not only stored but are also implicit in the position of the elements of the tree to be traversed.

Creation of the silicon substrate.
Octree construction is a well-known subject in the literature [10]. In our implementation, the only problematic issue is to dynamically generate the tree in combination with the use of a specific crystallographic unit cell that matches the desired orientation of the silicon crystal. To realize this objective, the crystallographic unit cell is created before generating the tree. This is done by creating a virtual crystal and finding potential sets of three orthogonal lattice vectors, with the only demand that they must satisfy the periodicity of the atomistic structure of the chosen surface, both in-plane and along the normal direction. The optimal orthorhombic unit cell is chosen by minimizing the aspect ratio of the two planar lattice vectors and minimizing the total cell volume. The atoms in the basis (i.e. inside the unit cell) are labeled without any special criterion and a list of the first neighbors for each atom in the basis is constructed. Using this orthorhombic unit cell, the tree can be constructed by using a standard recursive process. The orthorhombic unit cell is stored at the end of those branches that reach to the leaf level. The branches ending at lower levels represent empty regions, either in the bulk of silicon or in the etchant. For fast creation of every new tree node, we use bitwise concatenation of the tree coordinates.

Expansion and contraction.
Due to the fact that the root node of the tree defines the maximal spatial region for the simulation, the evolving surface may go out of the scope of the octree during the simulation. The number of cells (i.e. the tree leaves actually in use) can be used to repartition the space of the simulation, adding or reducing the number of levels in the tree. As an example, if n x , n y and n z are the dimensions along the x-, y-and z-coordinates of the tree and s = max(n x , n y , n z ), then the number of levels in the octree n is obtained as n = log 2 s. Note that it is not necessary that n x , n y and n z are equal, i.e. it is immaterial to make the space topologically cubic. When the surface gets dangerously close to the edge of the boundary, an expansion of the octree structure can be carried out to rebuild the system. If the tree expands once, the whole space can grow eight times larger than before. In this operation, the old root node simply converts to one child node in the new tree. On the contrary, if the whole bounding box around the system shrinks to only one-eighth of the defined space, the program can perform the inverse process of contraction to make computation more economical within a smaller scope.

Memory occupation.
The present CCA method benefits from the data representation in figure 2(a). During the simulation of the micro needle array shown in figure 1, the memory used by the octree structure is only 20% of the amount needed by the conventional dynamic method and the memory use shows an almost constant trend in the development of the 3D microstructure ( figure 2(b)).

Different types of material.
For the purpose of simulating the functionalities of different materials in MEMS applications, such as oxide, nitride and/or sacrificial layers, SOI wafers, etc, the octree can be enhanced by using different labels for atoms of different materials and controling the etching properties of the special atoms. The surface sites stored in the tree contain evolution information, including material properties, etch rates, as well as the position of the sites and the occupation state. This information is updated as the surface evolves.

Step flow
Prior to this study, typical methods such as KMC [12]- [14] and bond configuration model [15] have been applied to simulate anisotropic etching. Moreover different step flow models have been introduced to describe the anisotropy of the etch rate [16]- [18]. As figure 3 shows, in a simple step flow model the etch rate can be regarded as the composition of two processes with different velocities: propagation of the terraces and of the steps. This stems from the fact that the etchant attacks the step sites more actively than the terrace sites in such a way that a step atom is etched much faster than a terrace atom. If the difference between the rates of the two processes is large enough, etching occurs mainly at the edges of the steps. As a result, one may argue that the etch rate is primarily controlled by the density of steps.
In a previous model [16], the overall etch rate is geometrically expressed as a vector sum in terms of two velocities: where θ is the miscut angle from the close-packed plane. In the absence of a more detailed atomistic analysis, this kind of step flow model differs significantly from the experimental data, especially at certain orientations where the etch rate versus miscut angle has sharp discontinuities [16]. Additionally, this type of vector-based method is not directly applicable for the realization of atomistic simulations. However, it is possible to adapt the previous simplified step flow model for the atomistic scale and, in doing so, find a manner to implement the atomistic simulations. When we zoom into the atomistic structure of the steps of a generic vicinal orientation, usually four types of atoms can be distinguished, as shown in figure 3(b): the terrace atoms (terrace, T), the terrace atoms at the edge of the step (edge terrace, ET), the step atoms (step, S) and the terrace atoms restricted by the presence of the step (restricted terrace, RT). Typically, the terrace regions are always composed of (111), (110) and (100) pattern-like atoms. If the terrace width is large, the terrace contains examples of pure terrace sites, such as T for (775) in figure 3(b), while for shorter terraces there may not exist any pure sites, such as for (331) in figure 3(b), which only displays RT and ET sites.
In order to reproduce step flow, the ET and S atoms need to be very active while RT and T should be very stable. For instance, if S has the largest rate, step flow will proceed by removing S first and then ET (or, actually, a modified ET, since the type of the atom will change after removing S), repeating this periodically (i.e. S, mod-ET, S, mod-ET, S, mod-ET, etc). If ET has the largest rate, step flow is then realized as the removal of ET and modified S atoms (i.e. ET, mod-S, ET, mod-S, ET, mod-S, etc). Although the occupation of the sites on the terrace regions gradually decreases after each time step (section 2.2), eventually the removal of the corresponding atoms typically occurs only after they have been converted into one of the two active types on the step. This implies that the overall etch rate of the surface is essentially determined by the fractions of atoms on both the terrace and step regions, with the atoms on the step region, such as ET and S, playing the dominant role. In cases where the surface has wide terraces, the atoms on the terrace can be etched away before the propagating steps reach their positions. These terrace removals create other types of atoms, complicating the basic analysis.
Below, we will refer to the two active atoms at the steps as A and B, playing the roles of S and mod-ET, or ET and mod-S, depending on the specific realization of step flow.
In our model, we use a classification scheme for the surface atoms based on four indices (FS, FB, SS and SB), where FS and FB denote the numbers of first neighbors on the surface and in the bulk, respectively, and similarly SS and SB stand for the numbers of second neighbors on the surface and in the bulk, respectively [5]. The use of these indices allows to unambiguously characterize each and every of the surface sites. Although one can distinguish between more than 30 different sites by using the previous concepts of restriction (such as RT) and edging (such as ET) on vicinal (111), (110) and (100) together with the characterization scheme [9], it is not completely necessary to consider them all in order to understand the basic principle of the propagation of these surfaces according to step flow. More details about the definition of the sites are discussed in another paper [9].
The previous description using the removal of S and ET atoms is meant to be generic. When a particular vicinal orientation of (111), (110) and (100) is considered, certain specific sites play the roles of S and ET. Similarly, the crystallographic cut of a generic {hkl} surface can be considered as a collection of (111), (110) and (100) terraces or their combinations. See [19] for a physical insight on these aspects and [9] for a full detailed description of the different sites.

General considerations
The time evolution provided by the step flow CA model for the different crystallographic surfaces is characterized by the following features: (i) If the step propagation is dominated by the active atoms without any additional complication due to the removal of the terrace atoms, the etch rate is written as: where h j ( j = A, B) is the change in the average height of the surface when the atoms of type j on the step region are removed, and t j is the corresponding time increment. (ii) h j is a function of the miscut angle, which is different for each crystallographic plane.
See [9] for a full description of h j and the derivation of the values it takes for the different sites on the different surfaces. (iii) t j is related to the removal rate of the corresponding atom and the previous site-type history: where t is the total time spent by the atom as a previous type or types (typically as a restricted terrace and/or a terrace atom) and r j is the removal rate as an active atom on the step region (e.g. S, mod-ET, ET or mod-S). π j (t) = 1 − t o r (t )dt is the occupation of the site after exposure to the etchant for a time t. The term t o r (t )dt contains the past history of the atom before the removal. Using equation (3), the time increments of the two active 11 atoms satisfy the following equations: , ...} is the full site-type history of atom A, which is exposed to the etchant during several alternating time intervals of duration t A and t B . Thus, ) t B is the total reduction in the occupation of site A until atom A is removed. (The reduction in the occupation is necessarily 1.) In a similar manner, the second equation describes the total site-type history for atom B. For example, if the active step atoms are A = ET and B = mod-S, there are two equations for the time increments t ET and t mod−S . By solving equation 4, t ET and t mod−S can be expressed in terms of the rates for the atom types involved in the two site-type histories. See [9] for a detailed presentation of the particular site-type histories for a wide range of surface orientations and corresponding specific equations.
(iv) Once the height increments in (ii) and the time increments in (iii) have been expressed as functions of the geometrical measures and the removal rates of the involved atoms, the etch rate of each surface orientation in (i) can be formally expressed as a function of these geometrical and kinetic parameters. (v) In order to reduce the number of removal rates appearing in equation (4) and, eventually, in equation (2), we assume that the removal rate of the restricted atoms (such as RT) is equal to that of the non-restricted (such as T). Also, the removal rate of the modified-ET and S atoms (mod-ET and mod-S) is taken to be 1.0. This simplifies the final form of equation (2) for each surface orientation.

Time compensation (TC)
The TC method [3,5] has been introduced in order to improve the accuracy of the CCA simulations, which typically use constant time stepping (CTS). Since the occupation of a cell (or mass) is sequentially reduced in every time step by an amount equal to the (normalized) removal rate, the cell occupation eventually becomes smaller than the time step and only a segment of the time increment needs to be used when the cell is removed. TC adds the end part of the time increment to the next time step, thus compensating an otherwise erroneous measure of time. Unfortunately, TC works well only when the sites removed in a time step have an identical value of occupation and removal rate. TC does not correct the prevalent error associated with the use of CTS, namely, the simultaneous removal of different cells that should actually be removed sequentially according to their differences in occupation and rate. Thus, a simulation using TC does not necessarily follow the correct evolution history, e.g. for step flow, resulting in deviating etch rates for most orientations.
In this study, we use variable time stepping (VTS), always advancing time by an exact amount that matches the correct sequence, i.e. according to the values of the occupation and rate for all possible processes at the current state. This is very helpful for the derivation of equations (4) and (2) as it provides a manner to verify them. This exact time stepping is realized as follows. During a first loop over the whole surface, the minimum time to remove one surface 12 atom (or a group of atoms having the same time increment) is obtained as the time increment, so that only those atoms with this minimum value will be removed in a second loop over the surface. On the other hand, the use of carefully selected CTS (e.g. with a small time increment) together with TC, giving rise to TC-CTS, can provide similar output than exact VTS. The efficiency of a TC-CTS computation improves by using larger time steps but this may affect the accuracy and a trade-off is typically found. Both the exact (VTS) and the approximate (TC-CTS) methods are implemented in our program. VTS is used for deriving and verifying the analytical equations and TC-CTS is used for demonstrating the engineering applications.

Calculation of removal rates
It is possible to obtain values for the atomistic removal rates by inverting the relations obtained in section 3.2 (iv) for the etch rates of each crystallographic plane (R) and replacing the Rs by the experimental etch rate values. In this section, we describe a typical selection of planes that enables the construction of a system of equations, allowing to solve the atomistic removal rates. For presentation purposes, we focus on (111) (r 1 , r 2 , . . . , r n ), with i = A, B. Based on the sitetype definitions given in [9], the orientations in the range 0-54.7 • from (111) to (100) have five simplified step atom sites and two terrace sites. Thus, n = 7 in this case. In order to construct a system of equations for the seven unknowns, we need to use the etch rates of seven crystallographic planes, such as (111), (14 13 13), (533), (211), (311), (911) and (100). After that, the time increments can be substituted in equation (2), formally yielding expressions for the etch rate of each plane R {hkl} in terms of the removal rates {r 1 , r 2 , . . . , r n }. In this way, a set of seven non-linear equations is obtained, with the seven etch rates of the planes R {hkl} on the left-hand side and seven complex non-linear functions of the seven removal rates on the right-hand side. The atomistic removal rates r 1 , r 2 , . . . , r 7 are the unknowns.
Instead of solving the equations exactly, we seek the best fit to the experimental data by minimizing the sum of the squares of the residuals with respect to the experimental etch rates R exp {hkl} . We can state this optimization problem as: where r is the vector of unknowns r 1 , r 2 , . . . , r n , F(r ) = R exp − R cal (r ) is a vector and ε measures the global error. This optimization approach is necessary since the CA model is an idealization while many effects such as contamination, local etchant depletion and/or diffusion may result in deviations of the etch rates in the experiment.
After applying this procedure to a small set of orientations in the ranges (111)-(100), as described above, as well as (111)  It is important to stress the fact that our method provides the numerical solution for the etch rates by directly simulating the process on any crystallographic plane using the CA. This means that the CA can be directly used (without changes) to simulate the propagation of the etch front for engineering applications. This only requires adding the desired masking pattern and/or sacrificial layers to the desired wafer orientation, which is not limited to the principal surfaces. In comparison, previous step flow-based models devised for fitting the experimental etch rates are theoretical in nature, not allowing the direct simulation of the systems, especially for engineering purposes [16,17]. Moreover, the results of our study show that the exactsolution-assisted CA method is a flexible approach that can be applied to different etchants, at different temperatures and concentrations.

Stabilization of the etch rate
Because at the beginning, the occupation values of each surface atom is 1, the etch rate of any surface orientation shows an initial transient before it is stabilized. Since the active sites at the step dominate the surface propagation, the transient in the etch rate corresponds to the period required for the occupation field to reach a steady state over the surface. This behavior can be observed in figure 5, where a number of simulations were carried out for vicinal orientations of (111) differing in the miscut angle by increments of 3 • . The curves contain an initial exponential part and a stabilized constant part. The oscillation of the etch rate curves is due to the inherent periodicity of the step flow evolution provided by the CCA method. As considered above, the propagation of a step is split into the removal of two types of active atoms. The period of the oscillations corresponds to the time required for a step to propagate a distance equal to the width of the terraces, sometimes half of it, depending on the orientation. Because of this periodicity, the occupation field in the steady state is also periodic, not constant. The etch rates in figure 5 are averages, showing the oscillations only at the initial times. At later times the oscillations are averaged out.
The orientations with more active terrace sites, such as vicinal (110) and (100), typically present a more complicated evolution, not described by simple propagation of the steps. In these cases, the terrace regions are usually broken by the removal of fast terrace atoms. These faster etching planes create zigzag structures in the simulations. The (110) and (100) terrace-based orientations have the largest deviations from the simple step flow model. Independently of the simplicity or complexity of the evolution generated by the CCA method, the site-type history of any surface atom on any surface orientation experiences a well-defined sequence, similar to those considered in section 3.2 (iii).

Simulation results and discussion
The method presented in this paper significantly improves the description and simulation of the anisotropy of wet etching. The particular selection of orientations used in the computation of the non-linear least squares fit (equation (4)) may lead to different solutions for the removal rates in the system. The orientations should be selected under the consideration of minimizing the deviation from the experimental data. Figure 4 shows the fact that the anisotropy of the etch rate sometimes has a dramatic change as a function of the angle. The critical positions vary with different etchants and different temperatures. For example, a sharp discontinuity can be seen in the experimental curve at (551) (136.69 • ) in 50% KOH. A similar (although softer) discontinuity is observed at (533) (14.42 • ). Ignoring these points will produce a poor approximation.
Since the removal of the surface atoms is realized as a result of gradually decreasing the occupation of the sites, the presented step propagation model also allows the removal of pure terrace atoms, especially for those orientations presenting wide (110) or (100) terraces, where the terrace sites may be exposed for long times resulting in small occupation values. Typically, these removals occur close to the active step region, on the upper terrace. Even though this creates more complex morphologic features and evolution behavior, the surface can still reach a steady state. In particular, it can still be analyzed with the site-type history approach presented above, leading to non-linear equations and, correspondingly, it can be used finally in the optimization process. Instead of two equations for the two height and time increments of the two active atom types in the case of simple step flow, in these more complicated cases one can find three or more increments (in some cases even 8 or 10) with an equal number of equations. The selection of this type of surface will increase the difficulties in calculating and searching the optimized solution.
To demonstrate the practical interest of the developed simulation approach for engineering purposes, we consider the fabrication of a sharp micro needle array on a (100) substrate in 50% KOH at 70 • C, as shown in figure 6 [20]. This micromachining process is a combination of two dicing steps, realized by using a dicing machine, and one anisotropic wet etching step. In order to mimic the two dicing steps in the simulation, we simply etch vertically according to two mask patterns consisting of tiled squares. The process can be thought as two sequential deep reactive ion (DRI) etching steps. It generates square pillars with a smaller cross-section at the top due to the second step, as shown on the left-most frame of figure 6. After the vertical processing, maskless anisotropic etching is used to reshape the pillars into needle-like structures, as shown in the rest of the frames of figure 6. The description of the different stages of the process is very accurate in the simulations.
As a second example, figure 7 shows the micro fabrication of bridge-like underpass structures for micro fluidic channels [21]. The whole simulation involves the use of four  As a result of this study, a simulator of anisotropic and dry etching has been developed. Known as VisualTAPAS (Visual Three-dimensional Anisotropic Processing at All Scales), it can be freely downloaded from http://www.fyslab.hut.fi/∼mag/VisualTAPAS/Home.html. All the simulations presented in this study have been realized using this tool. Listings of the removal rates for different etching conditions can be found in the configuration file accompanying the program and can be easily accessed and/or modified (Simulation → Wet Etching → Rate settings → Advanced...). In addition to silicon substrates with virtually any orientation, one can also use SOI wafers, define sacrificial layers and both oxide and nitride masking patterns. The web site contains a comprehensive set of examples, including the details to reproduce the cases shown in figures 6 and 7.

Conclusions
This paper presents a new method for the simulation of complex fabrication processes involving anisotropic etching of silicon in MEMS. We show that, for this type of multivalued surfaces, the extended octree data structure not only represents efficiently the 3D topography of the dynamically evolving surface, but also shows better performance in the use of memory. The study shows that this advanced data structure can be well integrated into the underlying framework of a CA.
The combined use of the step flow model with a continuous description of the occupation of the surface sites using a CA enables a more accurate modeling of the anisotropic etching. This is partly a consequence of the derivation of exact, non-linear equations for the etch rates of the surfaces in terms of the atomistic removal rates. Moreover, a numerical solution is presented for fitting the atomistic rates based on solving the non-linear equations by minimizing the global error with respect to the experimental data. Based on a selection of some typical orientations, we can transform the macroscopic etch rates into values for the atomistic removal rates for direct use in the CA simulations. The positions where the etch rate experiences a sharp change as a function of the orientation play a crucial role in the determination of the atomistic removal rates.
The anisotropy of the etch rate can be well described and the simulated etch rates show good agreement with the experiments. Similarly, the simulation of complex MEMS applications is accurately performed. In addition, the simulation tool derived from the study can be used to gain a better understanding of the propagation of steps in anisotropic etching by showing how the different surfaces are etched as the etch front advances. The simulator has significant potential to model growth in micromachining also.