Comprehensive Simulations for Ultraviolet Lithography Process of Thick SU-8 Photoresist

Thick SU-8 photoresist has been a popular photoresist material to fabricate various mechanical, biological, and chemical devices for many years. The accuracy and precision of the ultraviolet (UV) lithography process of thick SU-8 depend on key parameters in the set-up, the material properties of the SU-8 resist, and the thickness of the resist structure. As feature sizes get smaller and pattern complexity increases, accurate control and efficient optimization of the lithography process are significantly expected. Numerical simulations can be employed to improve understanding and process design of the SU-8 lithography, thereby allowing rapid related product and process development. A typical comprehensive lithography of UV lithography of thick SU-8 includes aerial image simulation, exposure simulation, post-exposure bake (PEB) simulation, and development simulation, and this article presents an overview of the essential aspects in the comprehensive simulation. At first, models for the lithography process of the SU-8 are discussed. Then, main algorithms for etching surface evolvement, including the string, ray tracing, cellular automaton, and fast marching algorithms, are introduced and compared with each other in terms of performance. After that, some simulation results of the UV lithography process of the SU-8 are presented, demonstrating the promising potential and efficiency of the simulation technology. Finally, a prospect is discussed for some open questions in three-dimensional (3D) comprehensive simulation of the UV lithography of the SU-8.


Introduction
SU-8 is a negative, epoxy-type, ultraviolet (UV) photoresist based on EPON SU-8 epoxy resin from Shell Chemical, The Hague, The Netherlands, that has been originally developed, and patented (US Patent No. 4882245 (1989)), by International Business Machines Corporation (IBM, Armonk, NY, USA) [1]. The first SU-8 products were reported by MicroChem Inc., (previously Microlithography Chemical Corp., Westborough, MA, USA) [2]. Another company, Gersteltec Sarl, Pully, Switzerland [3], has also bought a license from IBM to produce and now sells the SU-8. Anyway, the products from the two companies are similar. Because thick SU-8 has excellent chemical stability, resistance against solvents, mechanical properties, and optical properties, it becomes a popular photoresist material for the fabrication of many microelectromechanical systems (MEMS) structures and devices [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18], including optical pick-up head [5], microlens [12], microfluidic channels [13], micropump [15], microneedles [17], cantilevers [18], and so on. Compared with X-ray LIGA (Lithografie, Galvanoformung, Abformung) [16] based on SU-8 using a very expensive synchrotron source [6,7], UV lithography of the SU-8 utilizes an inexpensive UV light source. This makes the UV lithography of the SU-8 more accessible. Currently, different UV lithography of SU-8, such as vertical lithography, inclined (tilted) lithography, multidirectional lithography, interference lithography, and so on, are widely used for fabricating MEMS structures and devices. It should be mentioned that there are several commercially available types of the same As shown in Figure 2, the aerial image simulation, exposure simulation, PEB simulation, and development simulation determine the final simulation profiles of the SU-8 for a typical comprehensive lithography simulation. For the beginning, the initial conditions and process parameters, including mask data, exposure time, UV light intensity, SU-8 layer thickness, UV light incident angle, gap between mask and SU-8 layer, reflectivity of substrate, PEB time and temperature for both two steps, development time and temperature, and so on, should be input into the simulation system. Then, the illumination of the mask on top of the SU-8 by the incident UV light source is obtained by aerial image simulation. During the following exposure simulation, the UV light propagation, as well as the generation of the photoacid HSbF6 within the SU-8, is simulated. After that, the post-exposure bake (PEB) simulation is used to describe the photoacid-catalyzed crosslinking reaction, obtaining the cross-linked site concentration in the SU-8 [56,57]. Thus, the development (etching) rate distribution into the SU-8 can be obtained according to the obtained cross-linked site concentration. Finally, with the development rate distribution, the development profiles of the SU-8 can be obtained by the appropriate etching surface evolvement algorithms introduced in Section 3. The final development profile data can be visualized and output for further analysis for various purposes. Generally speaking, the aerial image simulation and development simulation are time-consuming steps.  As shown in Figure 2, the aerial image simulation, exposure simulation, PEB simulation, and development simulation determine the final simulation profiles of the SU-8 for a typical comprehensive lithography simulation. For the beginning, the initial conditions and process parameters, including mask data, exposure time, UV light intensity, SU-8 layer thickness, UV light incident angle, gap between mask and SU-8 layer, reflectivity of substrate, PEB time and temperature for both two steps, development time and temperature, and so on, should be input into the simulation system. Then, the illumination of the mask on top of the SU-8 by the incident UV light source is obtained by aerial image simulation. During the following exposure simulation, the UV light propagation, as well as the generation of the photoacid HSbF 6 within the SU-8, is simulated. After that, the post-exposure bake (PEB) simulation is used to describe the photoacid-catalyzed crosslinking reaction, obtaining the cross-linked site concentration in the SU-8 [56,57]. Thus, the development (etching) rate distribution into the SU-8 can be obtained according to the obtained cross-linked site concentration. Finally, with the development rate distribution, the development profiles of the SU-8 can be obtained by the appropriate etching surface evolvement algorithms introduced in Section 3. The final development profile data can be visualized and output for further analysis for various purposes. Generally speaking, the aerial image simulation and development simulation are time-consuming steps. As shown in Figure 2, the aerial image simulation, exposure simulation, PEB simulation, and development simulation determine the final simulation profiles of the SU-8 for a typical comprehensive lithography simulation. For the beginning, the initial conditions and process parameters, including mask data, exposure time, UV light intensity, SU-8 layer thickness, UV light incident angle, gap between mask and SU-8 layer, reflectivity of substrate, PEB time and temperature for both two steps, development time and temperature, and so on, should be input into the simulation system. Then, the illumination of the mask on top of the SU-8 by the incident UV light source is obtained by aerial image simulation. During the following exposure simulation, the UV light propagation, as well as the generation of the photoacid HSbF6 within the SU-8, is simulated. After that, the post-exposure bake (PEB) simulation is used to describe the photoacid-catalyzed crosslinking reaction, obtaining the cross-linked site concentration in the SU-8 [56,57]. Thus, the development (etching) rate distribution into the SU-8 can be obtained according to the obtained cross-linked site concentration. Finally, with the development rate distribution, the development profiles of the SU-8 can be obtained by the appropriate etching surface evolvement algorithms introduced in Section 3. The final development profile data can be visualized and output for further analysis for various purposes. Generally speaking, the aerial image simulation and development simulation are time-consuming steps.

Aerial Image Simulation Models
Although the optical projection exposure method has been the main method for many years, contact and proximity exposure methods of the SU-8 using a so-called mask aligner are still the cost efficient exposure method in the current MEMS field [58,59]. The ultimate goal of aerial image simulation is to obtain the light intensity distribution in the SU-8. Because the wavelength of the 365 nm UV light is rather long, optical diffraction from masks will affect the precision of the lithography pattern. Furthermore, the reflection and refraction at the interface between the air gap (or compensation materials) and the SU-8, and the SU-8 and the substrate will also affect the precision of the lithography pattern [36,40,47]. Currently, aerial image simulations are based on two main categories of models. One category is based on scalar diffractive theory [29][30][31]38,41,47,59] and the other is based on rigorous electromagnetic field theory [60][61][62][63][64][65][66][67]. The former can be further divided into models based on the Fresnel diffraction [27][28][29][30][31]59], Fresnel-Kirchhoff diffraction [38,40,41,47], diffractive angular spectrum theory [33], and so on. While the latter can be further divided into models based on the finite difference time domain (FDTD) method [60][61][62] and waveguide method [63][64][65][66][67].
Compared with the models based on rigorous electromagnetic field theory [58][59][60][61][62][63][64][65], including the waveguide method, models based on the scalar diffractive theory are significantly faster (usually over hundreds of times or even more). Theoretically, these models are only valid for mask features that are large and thin compared with the wavelength of exposure light. This is not a problem for the UV lithography of the SU-8, which always satisfies this condition. The limitation of these models comes from the simulation accuracy, the application capability for special lithography methods such as inclined lithography, and the numerical costs for large-scale 3D simulations. To solve these problems, many studies have been conducted during the past several years. Tian et al. presented a correction to the Fresnel diffraction theory to increase the accuracy, depending on the scalar diffraction theory of light propagation [31]. For accurate and rapid simulation, incorporating the effect of reflection, transmission on the boundary in the oblique incidence case, and diffraction in the photoresist, Tang et al. extend the angular spectrum theory [68] to the light intensity distribution calculation for the lithography of thick photoresists. The propagation of the incident UV light will be refracted on the compensation material/SU-8 interface for the inclined UV lithography, and it is difficult to deal with the refraction in the Fresnel-Kirchhoff diffraction integral equation. To deal with this issue, Zhu et al. presented the "mask shifting approach" to handle the diffraction and refraction effects simultaneously [40]. To further reduce the light intensity distribution calculation time, Zhou et al. presented a method using an adaptable element size in x, y, and z directions for the whole calculation domain [47]. Using this method, the calculation time for 3D simulation can be reduced from 2 h to about 16 min for a moderately complex mask shape.
Compared with models based on scalar diffraction theory, the models based on rigorous electromagnetic field theory are more accurate. Since the pioneering work from Wong et al. [60], FDTD becomes an attractive method in rigorous mask diffraction simulation. FDTD uses a space domain discretization of the electric and magnetic field on equidistant staggered grids, solving the object fields at the mask rigorously. All diffraction, refraction, interference, absorption, and polarization effects are calculated out in the near field of the mask without approximation. FDTD can further reduce memory requirements and time per simulation by incorporating a graded mesh. The idea of the waveguide method is based on the theory of the Fourier modal method, and it is well adapted to lithography simulation with many typical mask geometries. The first prototype of the waveguide method was presented by Nyyssonen [63], and later this method was extended to a three-dimensional domain [64]. The concept of the waveguide method is the same. At first, the simulated structure is sliced into several layers with homogeneous optical properties along the vertical direction [29], z-direction, as shown in Figure 3. Then, the Fourier series is employed for the material parameters and electromagnetic field in all layers. Simultaneous equations formed by these expanded series and Maxwell Equations give an eigenvalue problem. After that, a hybrid approach using vector potentials is applied to couple electromagnetic fields in each layer. With extensive speed and convergence optimizations, an acceptable performance for this case is realized. Finally, using proper boundary conditions, the eigenvalue problem is solved, giving the light intensity distribution information throughout the propagation region.
Micromachines 2018, 9, x 5 of 21 using proper boundary conditions, the eigenvalue problem is solved, giving the light intensity distribution information throughout the propagation region. The combination of the waveguide method with domain decomposition techniques [66] and parallelized simulations enables the efficient and rigorous diffraction simulation for larger mask areas of the thin photoresist lithography [67]. The simulation time is less than 30 min on a 3.2 GHz Dual Core Pentium PC for a 245λ × 245λ × 6λ (λ is the wavelength of the exposure light) simulation areas. This means that a simulation speed over 15 times faster than that of the FDTD has been achieved. Theoretically, the models based on rigorous electromagnetic field theory can accurately simulate light intensity distribution in SU-8, if the limitations from the computation cost can be neglected. With the development of computer technology and the parallelization methods based on graphics processing units (GPUs), these limitations would seem to be partly solved in the near future.

Exposure Simulation Models
During the exposure process, the photoreaction initiator decomposes and a strong acid (HSbF6) is generated within the SU-8. According to the UV light intensity distribution into the SU-8, exposure models are used to describe the kinetics in the SU-8 during this process. The original Dill model [21], presented by Dill and his team at IBM for the mathematical equations describing basic lithography processes, has played an important role in conventional optical lithography simulations of thin photoresists for many years [68]. However, some special cases should be concerned for thick SU-8 compared with thin photoresists, such as the nonlinear factors in the exposure process of thick photoresists like SU-8. For example, the concentration distribution of contents and the refractive index vary along z direction, as shown in Figure 4. To solve this problem from the nonlinear characteristics, the original Dill model should be modified [33,42,69,70].  The combination of the waveguide method with domain decomposition techniques [66] and parallelized simulations enables the efficient and rigorous diffraction simulation for larger mask areas of the thin photoresist lithography [67]. The simulation time is less than 30 min on a 3.2 GHz Dual Core Pentium PC for a 245λ × 245λ × 6λ (λ is the wavelength of the exposure light) simulation areas. This means that a simulation speed over 15 times faster than that of the FDTD has been achieved. Theoretically, the models based on rigorous electromagnetic field theory can accurately simulate light intensity distribution in SU-8, if the limitations from the computation cost can be neglected. With the development of computer technology and the parallelization methods based on graphics processing units (GPUs), these limitations would seem to be partly solved in the near future.

Exposure Simulation Models
During the exposure process, the photoreaction initiator decomposes and a strong acid (HSbF 6 ) is generated within the SU-8. According to the UV light intensity distribution into the SU-8, exposure models are used to describe the kinetics in the SU-8 during this process. The original Dill model [21], presented by Dill and his team at IBM for the mathematical equations describing basic lithography processes, has played an important role in conventional optical lithography simulations of thin photoresists for many years [68]. However, some special cases should be concerned for thick SU-8 compared with thin photoresists, such as the nonlinear factors in the exposure process of thick photoresists like SU-8. For example, the concentration distribution of contents and the refractive index vary along z direction, as shown in Figure 4. To solve this problem from the nonlinear characteristics, the original Dill model should be modified [33,42,69,70].  The combination of the waveguide method with domain decomposition techniques [66] and parallelized simulations enables the efficient and rigorous diffraction simulation for larger mask areas of the thin photoresist lithography [67]. The simulation time is less than 30 min on a 3.2 GHz Dual Core Pentium PC for a 245λ × 245λ × 6λ (λ is the wavelength of the exposure light) simulation areas. This means that a simulation speed over 15 times faster than that of the FDTD has been achieved. Theoretically, the models based on rigorous electromagnetic field theory can accurately simulate light intensity distribution in SU-8, if the limitations from the computation cost can be neglected. With the development of computer technology and the parallelization methods based on graphics processing units (GPUs), these limitations would seem to be partly solved in the near future.

Exposure Simulation Models
During the exposure process, the photoreaction initiator decomposes and a strong acid (HSbF6) is generated within the SU-8. According to the UV light intensity distribution into the SU-8, exposure models are used to describe the kinetics in the SU-8 during this process. The original Dill model [21], presented by Dill and his team at IBM for the mathematical equations describing basic lithography processes, has played an important role in conventional optical lithography simulations of thin photoresists for many years [68]. However, some special cases should be concerned for thick SU-8 compared with thin photoresists, such as the nonlinear factors in the exposure process of thick photoresists like SU-8. For example, the concentration distribution of contents and the refractive index vary along z direction, as shown in Figure 4. To solve this problem from the nonlinear characteristics, the original Dill model should be modified [33,42,69,70].  Liu et al. presented an enhanced Dill model using a set of modeling parameters, taking into account of the effects caused by light diffraction or the scattering in the AZ4562 thick photoresist (MicroChemicals GmbH, Ulm, Germany) and the change of photoresist refractive index in the exposure process [71]. Similarly, Zhou et al. developed an improved Dill model describing the nonlinear effects in the SU-8 photoresist [45], to simulate the exposure kinetics during the exposure process by the following [38,47]: where the SU-8 was divided into n layers, and α i was the absorption coefficient for the ith layer. I 0 (x, y, z) was the light intensity. m(x, y, z, t) and C A (x, y, z, t) denoted the normalized photoacid generator concentration and HSbF 6 concentration, respectively. z(µm) was a variable to indicate the photoresist depth (the distance from top to bottom of the SU-8 layer). l was the distance that the UV light passes through in the SU-8. The Dill parameters varied with the SU-8 thickness: By fitting experimental results using the conventional "Dill graphical method" [72], the Dill parameters for various SU-8 thickness can be separately obtained, then i 1 , i 2 , i 3 , j 1 , j 2 , j 3 , l 1 , l 2 , and l 3 are fitted out. The Dill graphical method involves extracting the A and B parameters from the starting and ending sample transmittance values, while the C parameter is estimated from the slope of the initial portion of the transmittance curve. Based on this principle, some specific exposure parameter measurement tools have been available for many years, such as the ABC analyzer manufactured by Lithotech, Saitama, Japan [32,73,74]. For the SU-8 2000 series, the values for these constants were fitted to be as follows: i 1 = −0.0031, i 2 = 4.7878 × 10 −6 , i 3 = −3.0860 × 10 −9 , j 1 = 0.0079, j 2 = 7.5510 × 10 6 , j 3 = 4.0971 × 10 9 , l 1 = 0.0865, l 2 = −1.3078 × 10 −4 , and l 3 = 6.1193 × 10 −8 [38].

Post-Exposure Bake (PEB) Simulation Models
For the PEB process, HSbF 6 catalyzes ring-openings of epoxies and turns them into cross-linked sites. Thus, the SU-8 resin will transfer from a low-molecular-weight material to a highly cross-linked three-dimensional polymer network, and the degree of cross-linkage of SU-8 depends on the PEB temperature. The polymer network is significantly less soluble than the resin without cross-linking reaction. Bobbitt [75] introduces the details of catalytic chemical reactions of epoxy. Studies reveal that the physical properties, such as Young's modulus and the coefficient of thermal expansion, highly depend on UV exposure and PEB conditions [76], and these kinds of variations in physical properties have been investigated using molecular dynamics simulations in different groups [77,78]. However, this simulation approach can not currently be incorporated into the lithography profile simulations discussed in this paper.
The chemical reactions and diffusion of HSbF 6 happen simultaneously and couple with each other during the PEB process, because the HSbF 6 concentration is not uniform in SU-8. Thus, the coupled reaction-diffusion kinetics must be considered simultaneously. As for thin photoresists in the IC field, some complicated and accurate PEB models have been presented [24,25,62]. However, if they are expected to be used for the PEB process of the SU-8 lithography, many efforts should be put into the measurement of various model parameters. Currently, the model originally developed by Zuniga et al. [79] can be used to describe the PEB process for the SU-8. As the reaction order for HSbF 6 and cross-linked sites is 1 for the SU-8 lithography, the equations for the PEB process of the SU-8 can be expressed by the following [2,38,47]: where Equation (8) describes the reaction leading to the de-protection of the polymer sites, and Equation (9) describes the HSbF 6 evolution in time and space during the PEB process. C cs was the normalized concentration of the cross-linked sites where epoxy rings become open, and one cross-linked site is created when one epoxide ring is open. q was the order of the crosslinking reaction, D acid was the diffusivity of HSbF 6 , k 1 is the site-crosslinking reaction coefficient, and k 2 is the photoacid loss reaction coefficient. k 1 , k 2 , and D acid had an Arrhenius relationship with temperature. The coefficients in Equations (8) and (9)

Development Simulation Models
The first step to simulate the development of the SU-8 in developer is to calculate the dissolution (development) rate of each point in the SU-8. Various dissolution rate models have been proposed to describe the relationship between dissolution rate distribution within the SU-8 and cross-linked site concentration, such as the Mack model [82,83], Weiss model [84], Notch model [85,86], and Enhanced Notch model [85,86]. For example, the Notch Model is expressed as follows: where R min and R max are the minimum and maximum dissolution rate of the SU-8 photoresists, respectively; n 1 denotes the steepness of the rate curve; M th-notch is the threshold value of C cs where the notch occurs; and n n-notch is the measure for the strength of the notch. As shown by Equations (10) and (11), the parameter values for former mentioned models are dependent on the PEB conditions. The development parameters in the above-mentioned models can be determined using specific equipment, such as RDA-790 by Lithotech, Saitama, Japan [74], according to the reported procedures [30,87]. It should also be mentioned that some other parameters, besides the cross-linked site concentration, should also be measured for the above models. As we know, agitation methods are usually adopted for high aspect ratio structures during the development process, and the dissolution rate will be increased over several times, compared with the development process without an agitation method [2,38,47]. Furthermore, the dissolution rate of the SU-8 is depth-dependent [38,88], and the swelling effect will also affect the size of pattern edges and sidewall profiles of the SU-8 [38,42]. Although no accurate physical models for these effects are available, it is necessary to adopt the depth-dependent dissolution rate effect, and the swelling model, to improve the simulation accuracy for some cases [38,41,88]. We notice that the swelling effect seems more significant for the original SU-8 resist series than for its newer formulations (SU-8 2000 series and SU-8 3000 series) [2,44]. With the dissolution rate distribution in the SU-8, the propagation of the development profiles can be simulated using etching surface evolvement algorithms introduced in the Section 3.

String Algorithms
The string method was introduced by Jewett et al. [89] and was initially used to describe photoresist development for thin photoresists. The surface is represented by a string of nodes connected by straight line segments, as shown in Figure 5. The etching surface is advanced in small discrete time steps (∆t) by moving the etching point along a vector normal to the local surface. The advancement distance is the average of the advancement length of two adjacent segments. The string algorithms require little memory and computation time, as they only keep to track the etching surface. Jewett et al. [89] compared the string algorithm with the cell-based and ray-tracing algorithms using a test case of electron-beam exposure. They demonstrate that the string algorithm is more accurate and computationally more efficient than the cellular automata algorithms. The main disadvantage of the string algorithms is that loops will form in the boundary surface, especially when the etching rate in the material has a highly varied spatial distribution. The limited step size should be kept small to prevent loop formation and ensure accuracy of the solution. Segments must be removed when their length becomes comparable to the step size, otherwise loops can form. Loops form at the intersections of converging segments, which in turn results in physically incorrect surface evolution [94]. The de-looping algorithms must also be incorporated to ensure accuracy and correctness accuracy. Especially for the case of 3D simulations, the straight-line segments must be changed into triangles and polygons for representing the etching surface. Thus, the de-looping in 3D becomes much more computationally expensive, as it involves all triangle or polygon intersections in a mesh containing thousands of triangles or polygons.

String Algorithms
The string method was introduced by Jewett et al. [89] and was initially used to describe photoresist development for thin photoresists. The surface is represented by a string of nodes connected by straight line segments, as shown in Figure 5. The etching surface is advanced in small discrete time steps ( t Δ ) by moving the etching point along a vector normal to the local surface. The advancement distance is the average of the advancement length of two adjacent segments. The string algorithms require little memory and computation time, as they only keep to track the etching surface. Jewett et al. [89] compared the string algorithm with the cell-based and ray-tracing algorithms using a test case of electron-beam exposure. They demonstrate that the string algorithm is more accurate and computationally more efficient than the cellular automata algorithms. The main disadvantage of the string algorithms is that loops will form in the boundary surface, especially when the etching rate in the material has a highly varied spatial distribution. The limited step size should be kept small to prevent loop formation and ensure accuracy of the solution. Segments must be removed when their length becomes comparable to the step size, otherwise loops can form. Loops form at the intersections of converging segments, which in turn results in physically incorrect surface evolution [94]. The de-looping algorithms must also be incorporated to ensure accuracy and correctness accuracy. Especially for the case of 3D simulations, the straight-line segments must be changed into triangles and polygons for representing the etching surface. Thus, the de-looping in 3D becomes much more computationally expensive, as it involves all triangle or polygon intersections in a mesh containing thousands of triangles or polygons.

Ray Tracing Algorithms
A ray tracing algorithm was an application of the differential ray equation, and began to be used for etching surface advancement during lithography simulations in the 1970s [95]. In the ray tracing algorithms, the radial or directional vectors are initialized perpendicularly to the initial surface, and then each point moves along the direction vector in the discrete minute interval. After each step, the direction vector of the point is recalculated by the discrete ray equation. The developed and undeveloped regions are defined with different development rates. Actually, the propagation of the etching vectors is calculated like the propagation of light using Snell's law of refraction. The index of refraction of the SU-8 is defined as is the maximum value of the local development rate R(x, y, z). The evolvement of the etching front can be implemented based on the etching vectors. Because the ray trajectory is only determined by the previous trajectory and the corresponding etching rate at that point, there is no need to track the relationship between different rays. Theoretically speaking, the ray tracing algorithm is easy to implement. However, the actual situation is a little more complicated. The ray tracing algorithm does not describe the whole etching surface, but only the points on the surface. In order to reconstruct the etched morphology from the end point of the trajectory, a high density point is needed in a certain area. Even so, it is not easy to reconstruct the etching surface by the ray ends, because the rays are independent of each other.

Ray Tracing Algorithms
A ray tracing algorithm was an application of the differential ray equation, and began to be used for etching surface advancement during lithography simulations in the 1970s [95]. In the ray tracing algorithms, the radial or directional vectors are initialized perpendicularly to the initial surface, and then each point moves along the direction vector in the discrete minute interval. After each step, the direction vector of the point is recalculated by the discrete ray equation. The developed and undeveloped regions are defined with different development rates. Actually, the propagation of the etching vectors is calculated like the propagation of light using Snell's law of refraction. The index of refraction of the SU-8 is defined as n 1 = R max /R(x, y, z), where R max is the maximum value of the local development rate R(x, y, z). The evolvement of the etching front can be implemented based on the etching vectors.
Because the ray trajectory is only determined by the previous trajectory and the corresponding etching rate at that point, there is no need to track the relationship between different rays. Theoretically speaking, the ray tracing algorithm is easy to implement. However, the actual situation is a little more complicated. The ray tracing algorithm does not describe the whole etching surface, but only the points on the surface. In order to reconstruct the etched morphology from the end point of the trajectory, a high density point is needed in a certain area. Even so, it is not easy to reconstruct the etching surface by the ray ends, because the rays are independent of each other. Furthermore, the added rays must be in the sparse area of the ray, and sometimes some areas will not be etched because of inappropriate selection of the initial rays. The difficulties in reconstruction of the surface from independent rays make the ray tracing algorithms not very attractive for 3D etching simulations.

Cellular Automata Algorithms
The cellular automata algorithms were first presented by John von Neumann following Ulam's suggestions [98]. The computation domain is discretized into a set of square cells for 2D cases or cubic cells for 3D cases that are associated with different materials. The states at each cell are involved simultaneously based on the previous states of the cell and its neighboring cells at discrete time steps, according to a set of "local rules" [99,[102][103][104][105][106][107][108][109][110]. Because the "local rule" is only local relationships among neighboring cells, the governing equation for the whole computation domain is not necessary. Surface evolvement is represented by the removal or addition of cells or some parts of cells. For many years, the cellular automata algorithms have been used to simulate the shallow trench isolation etch with chlorine plasma [99,100], plasma etching of SiO 2 [101], anisotropic etching of silicon [101][102][103][104], photoresist development [105][106][107][108][109], and so on. A great advantage of the cellular automata algorithms is the ease of handling topological changes for arbitrary geometries and the easy with which they are implemented in 3D.
The preferential etching in different directions is great in the cellular automata algorithm with von Neumann neighborhood, significantly reducing the simulation accuracy compared with Moore neighborhood. For simplification, during the following part of this paper, discussions focus on the cellular automata algorithm with Moore neighborhood. For the development simulation, the SU-8 layer is divided into a matrix of identical square cells [106] or cubic cells [107,109] with side length a.
Here, 3D cases are discussed as an example. There are 6 adjacent cells, 12 diagonal cells, and 8 point cells in the neighborhood of a target cell. A cell is gradually etched by the etchant flowing from its neighbors. The local state of a cell is defined as the ratio of the etched volume to the total volume.
The 3D static cellular automaton algorithm is slow and inefficient [107], because the simulation program must process entire cells, even those lying in the interior of the photoresist, leading to reduction of the simulation speed and increase of the computer memory usage. For a simulation matrix of n × n × n, the computation time goes roughly as o(n 4 ), where n is the number of cells in one dimension. At the same time, the 3D static cellular automata algorithm needs to store the state array of the previous time step and the present time step (8n 3 or 16n 3 bytes). Furthermore, 4n 3 bytes are also required to store the etching rate array. Thus, the memory storage is larger than 12n 3 bytes. To improve simulation performance, Zhou et al. presented a 3D dynamic cellular automaton algorithm, using a dynamic memory allocation scheme [109], only processing the boundary cells. For the same simulation matrix of n × n × n, the computation time goes roughly as o(n 3 ), while the memory usage is about 6.4n 3 (2.4n 3 bytes to store the cell state flag array and the boundary cell array and additional 4n 3 bytes to store the etching rate array and the time compensation values). Detailed results, using a test function, revealed that for a simulation matrix of 100 × 100 × 100, the speed of the 3D dynamic cellular automaton algorithm is increased by over 10 times compared with that of the 3D static cellular automaton algorithm, while the memory usage is reduced by about 50% [109]. This makes the 3D dynamic cellular automaton algorithm practical for some 3D lithography simulation cases. Even so, the simulation speed is still a serious problem, limiting the application extension of the 3D lithography simulations of the SU-8. Here, a set of data is presented to illustrate this problem. For a simulation matrix with 500 × 500 × 500 using the 3D dynamical cellular automaton algorithm, the typical simulation time usually approaches several hours for a personal computer configuration: OS Windows XP SP3, CPU Intel Core2@2 GHz, DRAM 2 GB.

Fast Marching Algorithms
The fast marching algorithm, originally introduced by Sethian [111,112], briefly cited as fast marching method (FMM), is an optimally efficient algorithm for solving problems of front evolution where the front speed is monotonic. Namely, it is a stationary version of the general level set method [113,114]. The FMM is a powerful numerical technique for analyzing and computing moving fronts, and it has been applied to a wide variety of problems including seismic wave propagation, photoresist development, medical imaging, and optimal path planning. The flow of FMM can be divided into two parts: initialization and marching forward. In the initialization, the object is meshed into grids. During the marching forward to capture the etching surface evolvement, the FMM computes the time T (x, y, z) when the SU-8-developer interface passes through each point in a grid in the simulation area. The dissolution t rate of the SU-8, R (x,y,z) , is only a function of position, and the Eikonal equation ∇T (x,y,z) R (x,Y,Z) = 1 must be satisfied. Because the SU-8-developer interface can only advance in one direction during the whole photoresist dissolution progress, an "upwind" scheme for viscosity solutions of Hamilton Jacobi equations [115] can be adopted to solve the Eikonal equation.
Although the original FMM is accurate, fast, stable, and easy to implement in three dimensions, some research efforts have been made to further improve the original FMM by decreasing the time complexity and raising the method's precision. Yatziv et al. used a calendar untidy priority queue to reduce the time complexity [120]. Hassouna et al. adopted a multi-stencil to increase the precision [121]. From available results, the simulation speed of current fast marching algorithms was above eight times faster than the 3D dynamical cellular automaton algorithm [109]. However, original 3D FMM needs too many memory elements during simulations. For example, for a simulation matrix of n × n × n, the computation time goes roughly as 0(C*n 3 log(n)) + 0(n 2 ), where C is a pretty small constant and n is the number of grids in one dimension [109]. For memory elements, the original 3D FMM needs 12.4n 3 bytes of memory in total, including 4n 3 bytes to store the time values, 4n 3 bytes to store the dissolution rate array, 0.1 × 4n 3 bytes to store the min-heap structure of narrowband points, and 4n 3 bytes for companion array to point out the index of nodes in the min-heap. The simulation using the original FMM is faster than the 3D dynamical cellular automata algorithm, but needs too many memory elements. For the simulation matrix of 500 × 500 × 500 and the same personal computer configuration (OS Windows XP SP3, CPU Intel Core2@2 GHz, DRAM 2 GB) as in Section 3.3, the simulation cannot be implemented, as the memory limit of the computer is exceeded. To overcome this problem, Zhou et al. developed a 3D fast marching method based on full hash table (full hash fast marching method, briefly cited as HFMM) to reduce the memory usage [42].
As the original FMM only calculates nodes within the narrowband, a min-heap structure is developed to store the time values of the narrow band nodes. At the same time, a fast sweeping method is developed to search the node with the minimum time value in the narrow band, and to update the min-heap structure. Furthermore, a hash table is developed to store all data for the nodes in the narrow band, so the required memory elements are reduced. To illustrate the performance of the 3D HFMM, Figure 6 shows the running time and memory usage for a typical test function. Figure 6 indicates that the HFMM and the original FMM have nearly the same computational speed. However, nearly 60% of the computer memory elements in HFMM can be saved without any computational speed and accuracy loss, compared with the original FMM. With these advantages and the development of computer technology, algorithms extended from FMM will become increasingly competitive for the lithography simulation, especially for large-scale 3D lithography simulations of thick photoresists like SU-8.

Simulations and Discussions
Based on the models mentioned above and the algorithms for etching surface evolvement simulation, comprehensive simulation systems of the UV lithography process of SU-8 can be implemented. However, different level simulations of UV lithography, besides comprehensive simulations including all lithography steps or models, have been put into applications for various purposes and have achieved some success in many cases. For example, Cheng et al. [29] estimated the wall profiles and resolution for the near-field lithography of thick photoresist based on the Fresnel diffraction theory and exposure kinetics model. Chuang et al. [30] established an aerial image model to analyze the sidewall profiles of UV exposure of SU-8, leading to a novel method to reduce the diffraction problem from air gap using glycerol as a compensator to bridge out the air gap between the mask and SU-8. This method can greatly improve the sidewall straightness of high-aspect ratio SU-8 microstructures. Based on the effects of the Fresnel diffraction and absorption, Kang et al. [36] further presented an efficient model to predict the sidewall profiles of SU-8. Two exposure methods avoiding negatively sloped sidewall profiles were developed by utilizing the model, and SU-8 performs without undercut for the replication of optical waveguides were achieved using the exposure methods. Miao et al. [44] thoroughly analyzed and discussed the influences of different process parameters on the final surface profiles of micro lenses fabricated using UV lithography of SU-8, combining simulations and experiments. Yang et al. [46] built a complete 3D Fresnel-Kirchhoff diffraction model to predict the microstructure profiles fabricated by backside lithography of SU-8, combining a binary threshold approach. With the guidance of the above simulations, researchers can predict and optimize the lithography profiles [29,30], propose and evaluate some novel lithography techniques, and develop novel MEMS devices. However, it should be acknowledged that efficient comprehensive simulations are more suitable for understanding the fundamental effects of individual lithography step parameters and for optimizing the UV lithography process of SU-8. The reason lies in the fact that many lithography parameters will affect the final lithography profiles, including UV light distributions, exposure energy, development time, and so on. The accuracy of the simulation results is limited if very few parameters can be considered at the same time. For example, Figure 7 shows the simulation and experimental results for the lithography profiles of a cross-shaped mask with perpendicular incidence under UV source with 365 nm (2.6 mW/cm 2 ) radiation using the binary threshold approach. Here, the normalized light intensity of 0.33 is defined as the threshold value of SU-8. The maximum error between the simulation and experimental profiles in Figure 7 is about 10.4%. While the maximum error between the simulation profile of a comprehensive simulation process (Figure 8) and the corresponding experimental profile (Figure 7b) is about 7.6%. From this comparison, the disadvantage of the binary threshold approach in accuracy is evident.

Simulations and Discussions
Based on the models mentioned above and the algorithms for etching surface evolvement simulation, comprehensive simulation systems of the UV lithography process of SU-8 can be implemented. However, different level simulations of UV lithography, besides comprehensive simulations including all lithography steps or models, have been put into applications for various purposes and have achieved some success in many cases. For example, Cheng et al. [29] estimated the wall profiles and resolution for the near-field lithography of thick photoresist based on the Fresnel diffraction theory and exposure kinetics model. Chuang et al. [30] established an aerial image model to analyze the sidewall profiles of UV exposure of SU-8, leading to a novel method to reduce the diffraction problem from air gap using glycerol as a compensator to bridge out the air gap between the mask and SU-8. This method can greatly improve the sidewall straightness of high-aspect ratio SU-8 microstructures. Based on the effects of the Fresnel diffraction and absorption, Kang et al. [36] further presented an efficient model to predict the sidewall profiles of SU-8. Two exposure methods avoiding negatively sloped sidewall profiles were developed by utilizing the model, and SU-8 performs without undercut for the replication of optical waveguides were achieved using the exposure methods. Miao et al. [44] thoroughly analyzed and discussed the influences of different process parameters on the final surface profiles of micro lenses fabricated using UV lithography of SU-8, combining simulations and experiments. Yang et al. [46] built a complete 3D Fresnel-Kirchhoff diffraction model to predict the microstructure profiles fabricated by backside lithography of SU-8, combining a binary threshold approach. With the guidance of the above simulations, researchers can predict and optimize the lithography profiles [29,30], propose and evaluate some novel lithography techniques, and develop novel MEMS devices. However, it should be acknowledged that efficient comprehensive simulations are more suitable for understanding the fundamental effects of individual lithography step parameters and for optimizing the UV lithography process of SU-8. The reason lies in the fact that many lithography parameters will affect the final lithography profiles, including UV light distributions, exposure energy, development time, and so on. The accuracy of the simulation results is limited if very few parameters can be considered at the same time. For example, Figure 7 shows the simulation and experimental results for the lithography profiles of a cross-shaped mask with perpendicular incidence under UV source with 365 nm (2.6 mW/cm 2 ) radiation using the binary threshold approach. Here, the normalized light intensity of 0.33 is defined as the threshold value of SU-8. The maximum error between the simulation and experimental profiles in Figure 7 is about 10.4%. While the maximum error between the simulation profile of a comprehensive simulation process (Figure 8) and the corresponding experimental profile (Figure 7b) is about 7.6%. From this comparison, the disadvantage of the binary threshold approach in accuracy is evident.  Several research groups have reported their efforts to develop comprehensive simulation systems and to use the simulation systems for device design and process optimization. Zhou et al. have been developing comprehensive 2D and 3D simulation systems for UV lithography of the SU-8 for several years [37,41,44,47], and a series of 2D and 3D simulations have be implemented for vertical, inclined, and multi-directional lithography methods. Figure 9 shows the 2D simulation and experimental results of the inclined UV lithography of the SU-8 for 23.5° UV incident angle on bare silicon wafer [44]. Figure 9a,b reveal that the light absorption coefficient of the SU-8 is not a constant during the whole exposure processes. The reason lies in the fact that the cross-linked polymer has a higher absorption coefficient than the SU-8 before cross-linking reaction. The results also reveal that if only the UV light distribution is considered, the accuracy will be inevitably reduce. As both the exposed and unexposed SU-8 absorb material, the reflected UV light will be absorbed when passing through the SU-8. Thus, the energy of the reflected UV is relatively low to initiate cross-linking reactions of all SU-8, especially the SU-8 on the top, so the reflection effect is more notable at the bottoms of the SU-8 microstructures, as shown in Figure 9c,d. As more factors on the development profiles can be considered in the comprehensive simulation system [44], the simulation accuracy has been increased, compared with the profiles obtained using aerial image simulation and the threshold approach [40]. The largest line width variation between the experimental and simulation profiles in Figure 9c,d is less than 1.72 μm. Figure 10 shows the 3D simulation and experimental  Several research groups have reported their efforts to develop comprehensive simulation systems and to use the simulation systems for device design and process optimization. Zhou et al. have been developing comprehensive 2D and 3D simulation systems for UV lithography of the SU-8 for several years [37,41,44,47], and a series of 2D and 3D simulations have be implemented for vertical, inclined, and multi-directional lithography methods. Figure 9 shows the 2D simulation and experimental results of the inclined UV lithography of the SU-8 for 23.5° UV incident angle on bare silicon wafer [44]. Figure 9a,b reveal that the light absorption coefficient of the SU-8 is not a constant during the whole exposure processes. The reason lies in the fact that the cross-linked polymer has a higher absorption coefficient than the SU-8 before cross-linking reaction. The results also reveal that if only the UV light distribution is considered, the accuracy will be inevitably reduce. As both the exposed and unexposed SU-8 absorb material, the reflected UV light will be absorbed when passing through the SU-8. Thus, the energy of the reflected UV is relatively low to initiate cross-linking reactions of all SU-8, especially the SU-8 on the top, so the reflection effect is more notable at the bottoms of the SU-8 microstructures, as shown in Figure 9c,d. As more factors on the development profiles can be considered in the comprehensive simulation system [44], the simulation accuracy has been increased, compared with the profiles obtained using aerial image simulation and the threshold approach [40]. The largest line width variation between the experimental and simulation profiles in Figure 9c,d is less than 1.72 μm. Figure 10 shows the 3D simulation and experimental Several research groups have reported their efforts to develop comprehensive simulation systems and to use the simulation systems for device design and process optimization. Zhou et al. have been developing comprehensive 2D and 3D simulation systems for UV lithography of the SU-8 for several years [37,41,44,47], and a series of 2D and 3D simulations have be implemented for vertical, inclined, and multi-directional lithography methods. Figure 9 shows the 2D simulation and experimental results of the inclined UV lithography of the SU-8 for 23.5 • UV incident angle on bare silicon wafer [44]. Figure 9a,b reveal that the light absorption coefficient of the SU-8 is not a constant during the whole exposure processes. The reason lies in the fact that the cross-linked polymer has a higher absorption coefficient than the SU-8 before cross-linking reaction. The results also reveal that if only the UV light distribution is considered, the accuracy will be inevitably reduce. As both the exposed and unexposed SU-8 absorb material, the reflected UV light will be absorbed when passing through the SU-8. Thus, the energy of the reflected UV is relatively low to initiate cross-linking reactions of all SU-8, especially the SU-8 on the top, so the reflection effect is more notable at the bottoms of the SU-8 microstructures, as shown in Figure 9c,d. As more factors on the development profiles can be considered in the comprehensive simulation system [44], the simulation accuracy has been increased, compared with the profiles obtained using aerial image simulation and the threshold approach [40]. The largest line width variation between the experimental and simulation profiles in Figure 9c,d is less than 1.72 µm. Figure 10 shows the 3D simulation and experimental profiles of inclined UV lithography for 30.0 • incident angle, with TiO 2 film as an antireflection layer. The simulation time is about 39 min for the whole lithography process (an impressive high simulation speed), because adaptable element size and HFMM are adopted [47]. The maximum error between the simulation and experimental profiles in Figure 10 is less than 7.3%. With the simulation system, Zhou et al. have investigated the effects of different lithography parameters on the line width and topography deviation for lots of SU-8 microstructures, including exposure time, development time, line/space, and mask shapes [39,47].
Compared with current lithography simulation tools of thin photoresists in the IC field, the applications of UV lithography simulations of the SU-8 or other thick photoresists are significantly limited. The lithography simulations are now critical and standard tools for IC fabrication, thin photoresist design, and its evaluation. To some degree, UV lithography simulations of the SU-8 are behind of the lithography technology development of the SU-8. Up to now, few simulation tools are available. IntelliSense, a famous MEMS CAD company (Nanjing, China), released a simulation module Exposure for the UV lithography of SU-8 by collection with Southeast University, China [122].
profiles of inclined UV lithography for 30.0° incident angle, with TiO2 film as an antireflection layer. The simulation time is about 39 min for the whole lithography process (an impressive high simulation speed), because adaptable element size and HFMM are adopted [47]. The maximum error between the simulation and experimental profiles in Figure 10 is less than 7.3%. With the simulation system, Zhou et al. have investigated the effects of different lithography parameters on the line width and topography deviation for lots of SU-8 microstructures, including exposure time, development time, line/space, and mask shapes [39,47].
Compared with current lithography simulation tools of thin photoresists in the IC field, the applications of UV lithography simulations of the SU-8 or other thick photoresists are significantly limited. The lithography simulations are now critical and standard tools for IC fabrication, thin photoresist design, and its evaluation. To some degree, UV lithography simulations of the SU-8 are behind of the lithography technology development of the SU-8. Up to now, few simulation tools are available. IntelliSense, a famous MEMS CAD company (Nanjing, China), released a simulation module Exposure for the UV lithography of SU-8 by collection with Southeast University, China [122].  [44].
To extend the application of lithography simulations of the SU-8, several issues should be concerned. From the view of practical application, 3D simulations for various masks, as shown in Figure 11, even full chip level simulations, are greatly expected for the lithography process optimization and dimensional control of microstructures. However, these simulations require too many computer resources, and the computer resources for certain PCs are certainly limited. To deal with this problem, fast and efficient algorithms should be developed. Furthermore, current main stream algorithms, such as cellular automaton and fast marching algorithms, have an intrinsic parallel updating nature, and the corresponding simulations are highly inefficient when performed To extend the application of lithography simulations of the SU-8, several issues should be concerned. From the view of practical application, 3D simulations for various masks, as shown in Figure 11, even full chip level simulations, are greatly expected for the lithography process optimization and dimensional control of microstructures. However, these simulations require too many computer resources, and the computer resources for certain PCs are certainly limited. To deal with this problem, fast and efficient algorithms should be developed. Furthermore, current main stream algorithms, such as cellular automaton and fast marching algorithms, have an intrinsic parallel updating nature, and the corresponding simulations are highly inefficient when performed on classical central processing units (CPUs). Thus, improved parallel methods for the UV lithography simulation of the SU-8 on graphics processing units (GPUs) will be helpful. The recent cost-efficient examples of the parallel architectures for silicon anisotropic etching simulations based on GPU [123,124] provide a useful reference for this kind of research. Similarly, former mentioned image simulations based on rigorous electromagnetic field are highly flexible and accurate, and deal relatively easily with various problems in lithography simulation. However, the computer calculation ability limits their applications in UV lithography simulation of the SU-8. The usage of GPUs is also a possible solution to this problem using massively parallel computing approaches. In addition, an important aspect for a lithography simulator is accuracy. More accurate (if possible) physically-based models for various lithography simulation steps are expected over a wider range of process conditions. Finally, more efforts are expected to extract and optimize more model parameters. Research efforts along these lines will make the simulation technology more practical and this can satisfy the needs of UV lithography technology development of the SU-8.  [123,124] provide a useful reference for this kind of research. Similarly, former mentioned image simulations based on rigorous electromagnetic field are highly flexible and accurate, and deal relatively easily with various problems in lithography simulation. However, the computer calculation ability limits their applications in UV lithography simulation of the SU-8. The usage of GPUs is also a possible solution to this problem using massively parallel computing approaches. In addition, an important aspect for a lithography simulator is accuracy. More accurate (if possible) physically-based models for various lithography simulation steps are expected over a wider range of process conditions. Finally, more efforts are expected to extract and optimize more model parameters. Research efforts along these lines will make the simulation technology more practical and this can satisfy the needs of UV lithography technology development of the SU-8.    [123,124] provide a useful reference for this kind of research. Similarly, former mentioned image simulations based on rigorous electromagnetic field are highly flexible and accurate, and deal relatively easily with various problems in lithography simulation. However, the computer calculation ability limits their applications in UV lithography simulation of the SU-8. The usage of GPUs is also a possible solution to this problem using massively parallel computing approaches. In addition, an important aspect for a lithography simulator is accuracy. More accurate (if possible) physically-based models for various lithography simulation steps are expected over a wider range of process conditions. Finally, more efforts are expected to extract and optimize more model parameters. Research efforts along these lines will make the simulation technology more practical and this can satisfy the needs of UV lithography technology development of the SU-8.

Conclusions
The essential aspects in the comprehensive simulation of UV lithography of the SU-8 are reviewed and discussed. This shows that great progress has been made in this field during the past several years. Available lithography models for exposure, PEB, and development steps are analyzed, and the main algorithms for the etching surface evolvement simulation are reviewed with their advantages and challenges. Some typical simulations of the UV lithography of the SU-8 are also introduced and discussed. Finally, several challenging issues regarding the simulation speed and accuracy are discussed, to realize further extension of the applications of the UV lithography simulations of the SU-8.