{'fraction_non_alphanumeric': 0.032315013816395455, 'fraction_numerical': 0.008443352778630642, 'mean_word_length': 5.040333796940195, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical simulation of surface diffusion motion and its application in MEMS fabrication

Surface diffusion motion has been proposed as a new technique, called Silicon Migration Technology (SiMiT) for the fabrication of Micro-Electro-Mechanical System (MEMS) devices. Using this new technology, the integration of MEMS-CMOS can be achieved without the “sacrificial-structural” layers and “release-etch” procedure, which are necessary in traditional approaches. However, numerical simulation of surface diffusion process is still a challenging task due to the stiffness nature and the topological changes during the process. In this paper, we adopted a high-order explicit level-set method based on the total variation diminishing Runge-Kutta (TVD-RK) method, a high-order scheme for distance computation. An efficient smoothing scheme is proposed to increase the numerical stability and an adaptive time marching scheme is developed to accelerate the simulation. Experimental investigation of shape transformation of silicon trenches is carried out to verify the accuracy of the numerical scheme. Results show that numerical simulation results match quite well with experimental ones, which suggests the accuracy and reliability of the numerical simulation in the application of MEMS design.


Introduction
MEMS devices have been widely used in a variety of fields, including, but not limited to, fluidics, robotics control, biomedical industry and optical displays [1][2][3][4][5][6][7][8][9][10][11]. However, monolithic integration of MEMS with driving, controlling, and signal processing electronics is necessary in the manufacture procedure to achieve miniaturization and improve the performance, which is still a challenging task. Traditional MEMS-CMOS integration approach can be categorized into pre-CMOS, intra-CMOS and post-CMOS [12; 13]. In pre-CMOS process, the MEMS structures are released by a wet etch after the completion of the standard CMOS processes. This process is made complicated by the need of a planarization step and the thickness of the MEMS structures is constrained by the lithographical limit. The dimension limit also exists in the intra-CMOS process. High cost and the need for dedicated manufacturing foundries limit the application of pre-CMOS and intra-CMOS [12]. In the post-CMOS process, careful selection of a release-etchant with good selectivity between the structural and sacrificial layers is required [14; 15]. The necessary of "sacrificial-structural" layers and "release-etch" procedure for the traditional integration approaches increases process complexity and limits device packing density. To solve this issue, recently, a new technique based on surface diffusion of silicon is proposed and utilized in the integration of MEMS-CMOS without a sacrificial layer etch [16][17][18]. A buried cavity sealed with a silicon cover-diaphragm is first formed. The electronic devices, both on and off the diaphragm, are next fabricated. With a sealed cavity, the issues of material-and processincompatibility inherently present in earlier integration schemes are largely avoided.
This new technique, however, has not been widely used because of a lack of fundamental understanding of the phenomenon and a lack of design tools. To understand the phenomenon, predict the evolution of the micro structures and guide the experiments, we developed an efficient and accurate numerical method based on level-set method to investigate the silicon migration process. The level-set method was first proposed by Osher and Sethian [19] to represent interfaces implicitly. The most significant advantage of the level-set method lies in its ability to capture topological changes, such as merging or breaking, in a natural way. This method has been widely applied in a variety of areas including, but not limited to, multi-fluid flows [20][21][22], image process [23][24][25], film and crystal growth [26; 27]. Recently, it has been extended to modeling surface diffusion motion [28][29][30][31].
Direct simulation of surface diffusion is tough due the stiffness nature of the governing equation. In [28], Smereka proposed a semi-implicit scheme in which the stiff linear part is extracted and treated implicitly, while the nonlinear part is treated explicitly. This method was further enriched by adding the gradient information to obtain the subgrid resolution [29]. In [30], Bruchon proposed a mixed formylation to stabilize the numerical scheme. However, for these methods, the temporal accuracy is only first-order. The extension to higher-order scheme is yet to be developed. In addition, reinitialization and extension procedures are necessary, which may cause error accumulation, when solved using standard PDE-based method.
In this paper, a high-order level-set method is proposed for the simulation of surface diffusion motion. Experiments are also conducted by investigating silicon annealing under high temperature. Good agreement with experimental results suggests the accuracy and reliability of the numerical simulation in the application of MEMS design.

Methodologies
In this paper, the high-order level-set method proposed in [31] is adopted for the simulation of the surface diffusion process. For efficient consideration, a smooth procedure different from [31] is adopted and an adaptive time scheme is also utilized to further increase the efficiency. Finally, the experimental setup is presented.

Physical model
Moving interfaces can be described in the following form by the Lagrangian representation (1) where represents the evolving surface; is the normal velocity; is the unit normal vector. The normal velocity via surface diffusion is given by [32] ( where is the mean curvature; is a material constant that is proportional to the surface diffusivity and surface tension; , the surface laplacian, is the so-called Laplace-Beltrami operator. Eqs. (1) and (2) can be solved directly by front tracking method [33], through discretizing the interfaces, calculating the normal velocity and evolving the interfaces. However, this method suffers from the difficulty in handling complex topological changes such as pinch-off and merging. In this paper, an implicit-interface method (level-set method) is used instead to solve the surface diffusion problem.

High-order level-set method
The level-set method was proposed by Osher and Sethian [19] as an implicit fronting tracking method and recently, it has been extended to surface diffusion simulation [28][29][30]. In order to improve the accuracy and capture complex topologic changes, a high-order level-set method is developed in [31] to perform reinitialization and extension. The high-order scheme basically consists of four steps (as shown in Figure 1  As the velocity is only defined on the interface, for any point in the simulation domain, the velocity is set as the velocity of its closest point on the interface, which is calculated using the polynomial representation of the level-set function within the cell containing . In this way, the velocity extension is performed automatically, as the velocity is evaluated at the closest point on the interface.

Smoothing procedure
For surface diffusion motion, the velocity involves forth derivatives of the level-set function, which is prone to the numerical error. In addition, the high-order closest point algorithm is essentially entirely geometric, which results in a noisy velocity field. To ensure the stability of the numerical scheme, an accurate and smooth velocity field is required. In this work, the following regularization is adopted for the smoothing of curvature and normal velocity: where and are the curvature and normal velocity calculated by the method discussed previously, and are the smooth curvature and normal velocity after the regularization, is the time step and and are positive constants. In this work, a constant value 0.5 for both and is used for all the simulations. Finally, the level-set function will be updated by third order total variation diminishing Runge-Kutta scheme [34]. This scheme can guarantee that no spurious oscillations are produced as a consequence of the higher order accurate temporal discretization as long as no spurious oscillations are produced with the forward Euler building block [35].

Adaptive time step
In this work, instead of using a constant time step, we developed an adaptive time scheme based on the maximum of the normal velocity. The basic principle of the choice of the adaptive time is to guarantee that within one-time step, the interfaces will not move more than one grid cell. This choice will prevent large time steps when the interfaces change significantly and allow large time steps when the change of the interfaces is smooth, which can help to capture the interface accurately and efficiently. This kind of adaptive scheme is especially useful for capturing topological changes. For example, when there is no topological change, the interface will usually change smoothly and can be captured using a larger time step. When topological change happens (merging or self-splitting), sharp corners will appear, which renders significantly large velocity, and a smaller time step is required.
To determine the time step, we first obtain the maximum of the normal velocity and then calculate . The time step is set to be where is used to control the time step, and are used to prevent the time step from being too small or too large. In our simulations, , and . All these quantities are in non-dimensioned form and the corresponding equation is given by Eq. (6) and you can also refer to [31] for more information. Here we should notice that even though this kind of choice violates the CFL restriction ( ), the numerical is still stable, because high-frequency oscillations can be filtered out by the regularization (Eqs. (3) and (4)).

Experimental setup
Experimental investigation of shape transformation of silicon trenches during annealing is also studied. Conventional p-type silicon wafer is used in this study. An array of trenches (Figure 2 (a)) is first fabricated using a deep reactive-ion etcher (DRIE). Then heat-treatment is carried out at in argon at the atmospheric pressure to form the final cavity (Figure 2 (c)). By varying the initial configuration (parameterized by the array period , the diameter and the depth ), final cavities with different sizes and locations can be obtained. Interested readers can refer to [17; 18; 36] for more information.

Results and discussion
In this part, we studied temporal evolution of 3D rectangular trench arrays both numerically and experimentally. The initial configuration of the trench is illustrated in Figure 3. , and are set to be , and separately. The comparison of trench profiles between the numerical and experimental results at different time is shown in Figure 4.   It is shown that if the bottom part of the interface is cut from an inclined plane with from the vertical plane, the numerical results match quite well with the experimental results. This is actually due to the fact that in the experiment, it is very hard to guarantee that the cutting plane of the cross-section is vertical, as is shown in Figure 5. It is obvious that in the bottom part, the cutting plane is not vertical.
For numerical simulation, we choose to be one without loss of generalization. From dimensional analysis, we can obtain that ⁄ The obtained parameter is the same order with the estimated value from physical parameters in [33].
The efficiency of the adaptive time step is also studied. Figure 6 shows the time step as a function of the evolution time in the simulation of periodic 3D trench array. Initially, the time step is very small because of the corners existing in the initial trench array. Then the time step becomes larger and larger as the interface becomes smoother and smoother until the self-splitting topological change happens at the top of the trench. Before the splitting moment, the time step drops and after the splitting, the time step increases gradually. Before the second topological change (merging) happens, the interface is smooth and thus, the time step is larger. At the moment of merging, sharp corners are formed near the merging points and thus, the time step is dropped suddenly. After the merging, the time step begins to increase again. Thus, the adaptive time scheme can adjust the time step according to the smoothness of the interface and thus the simulation can be performed more efficiently without significant loss of accuracy.

Summary
In this paper, a high-order level-set method is adopted for the simulation of surface diffusion. A smoothing procedure and adaptive time marching scheme are proposed to increase the numerical efficiency. The shape evolution of silicon trenches during annealing is investigated both numerically and experimentally. Good agreement with experimental results also suggests the accuracy and reliability of the numerical simulation in the application of MEMS design.