1. Introduction

Laser direct writing lithography (LDWL) has been also widely used in the fields of 3D micro-optical elements [1,2], cell cultures scaffolding [3], micro-electro-mechanical system [4,5], and microfluidic chips in recent years as a maskless lithography technology to counter the increase of semiconductor industry mask costs, owing to the high degree freedom in controlling the element morphology [6,7]. However, manufacturing challenges still exist, such as the optical proximity effect (OPE) [8], caused by the LDWL system’s operation and the more significant optical absorption phenomenon of a thick-film photoresist compared with a thin-film photoresist [9], resulting in the morphology-fidelity reduction of manufactured 3D micro structures.

The optical proximity effect correction (OPC) method is well-known for improving image fidelity in integrated circuits. The thin-film photoresist model is used in the OPC process, and the shape and placement errors of 2D manufacturing targets are corrected. However, in contrast to the OPC method, which aims to control the line width and position accuracy of semiconductor devices, the 3D OPC method aims to control the 3D morphology of micron-scale structures. Because morphology is important in some cases, for example, aspheric lenses have a better aberration reduction effect than spherical lenses [10], and some free-form surface devices will also be designed to achieve the effect of beam shaping [11]. However, the realization of these functions depends on the morphology of free-form surface designed by researchers. Therefore, making microstructures with free-form surfaces is an important challenge to realize designed functions. In this study, to achieve better control of 3D micro device morphology in laser direct writing, a 3D OPC method is proposed, considering both the OPE and nonlinear response of thick-film photoresist (TF-NLR).

To overcome the OPE, it is necessary to recalculate the exposure dose distribution on the upper surface of the photoresist. Because of LDWL’s point-by-point exposure point-by-point optimization is possible. The aerial energy calculation process of LDWL is the superposition of the illumination intensity of adjacent exposure positions under the assumption of incoherent system, which is similar to convolution and will be discussed in the process of modeling 3D LDWL.

Another key step is to identify the characteristics of TF-NLR. Actually, predecessors have done extensive research in solving and simulating the properties of photoresist using two methods: layered numerical models and exposure-development curve fitting [12,13]. For the former, considering the change of refractive index of photoresist during exposure, energy absorption in photoresist, and concentration distribution of photosensitizer, the developed morphology can be predicted theoretically [14–16], because the local development speed determined by the concentration distribution of photosensitizer can be integrated with time to get the local development depth. This ensures that more parameters, such as developing time, can be optimized in the optimization process, but iteration and numerical simulation increase the computational complexity. To reduce the computational cost, we adopt the other method [17]. Referring to the development time when the side wall angle is the largest in research [18], by collecting the nonlinear relationship between the aerial energy on the photoresist surface and the local development depth of the photoresist in a pre-experiment, the expected morphology can be obtained by inverse optimization of the aerial energy with low computational cost. The study shows that 3D OPC method has achieved a good effect on the control of morphology under the recommended developing process. This method can be extended in the future. By comparing the optimization effects of different development times, the development time can be optimized in subsequent study.

Currently, in the process of dose distribution optimization, interior point method [19], steepest descent method [20], constrained gradient search method (L-BFGS) [21], Lagrange multiplier method, and deep learning have been applied [22,23]. In this case, when the step size can be adjusted automatically, local minima should be avoided. After comprehensive consideration, Adam gradient descent method, which is commonly used in deep learning training, was adopted as the optimizer [24]. The result showed that it could solve the dose distribution effectively. Hereto, the 3D OPC method based on LDWL has been established. Combining the 3D LDWL model with the properties of thick-film photoresist, a 3D target structure is input into the optimizer, and the optimized numerical dose distribution is provided after iteration in optimizer. The framework is shown in Fig. 1.

 

Fig. 1. 3D OPC method optimization framework.

Download Full Size | PDF

To verify the effectiveness of 3D OPC method in the manufacture of 3D free-form morphology structures and stability of the optimizer, we conducted an experiment on a group of quadratic curvature microlens arrays as an example. The result and analysis of the experiment show that the structures morphology has been effectively controlled.

2. 3D OPC method

2.1 3D LDWL model

When establishing the physical model of the 3D LDWL system, the single laser beam emitted by the light source was simplified into a Gaussian beam propagating along the Z direction (Fig. 2(a)), in which the waist plane of the beam (Z = 0) is placed at a fixed position relative to the photoresist surface, and the actually measured light intensity distribution on the waist plane of the beam in Fig. 2(a) was fitted as a Gaussian distribution (Fig. 2(b)) which is expressed by Eq. (1). Some assumptions are required. The focal plane of the laser beam in Fig. 2(a) is kept stably focused on the upper plane of the photoresist, and structures are fabricated within the depth of focus (DOF) range. For the 3D LDWL model, under the assumption of a linear time-invariant system (LTI-system) and incoherent imaging, in the case of two-dimensional distribution, the aerial image distribution $\textrm{E}(x,\;y)$ is obtained by convolution of light source ${G(x,\;y)}$ and dose distribution ${D(x,\;y)}$ in Eq. (2), and the aerial energy $\textrm{E}(x,\;y)$ is mapped by photoresist function (p function) to obtain the developed morphology ${M(x,\;y)}$ expressed in Eq. (3). The above flow is shown in Fig. 2(c).

 

Fig. 2. a) Sliced 3D model of Gaussian beam, which is simplified from light source beam, the coordinate axis unit is 100 nm/pixel. To simplify calculation, the light intensity distribution measured at the waist of the beam is taken as a 2D light source model and fitted by a Gaussian distribution. b) Normalized light source intensity distribution model based on the Gaussian distribution hypothesis (2D and 3D), which is the 850 nm (full width at half maxima, FWHM) laser spot data obtained by actually measuring the light source of our equipment. It has been processed by Gaussian smoothing, and the grid is divided into 100 nm/pixel. c) Flow chart of LDWL based on the thick-film photoresist. As an example, the fabrication process of microlens with 3D morphological requirements is provided. Developed morphology ${M}({x,y} )$. will be different based on whether the photoresist is positive or negative.

Download Full Size | PDF

Here, ${(x,\;y)}$ is the position coordinate of each point in the working space, ${w_0}$ denotes the spot size, and a denotes the coefficient scalar of the light source model. And $\otimes $ is the convolution operator.

2.2 Thick-film photoresist characteristic

We assume that the thickness of the photoresist is uniform. Furthermore, the reflection of the photoresist surface and substrate as well as the defects of photoresist surface are neglected. Because of the nonlinear factors, such as the absorption of thick-film photoresist, dynamic changes in refractive index during exposure, and uneven distribution of photosensitive agent concentration in the thickness direction caused by pre-baking [25], the output energy of laser does not linearly correspond to the actual development depth. However, when compared to the LDWL convolution model, thick-film photoresist plays a significant role in the manufacture of 3D structures, and it is also the determining factor to whether theoretical dose optimization can be applied to practice. The issue is that the numerical solution of photoresist properties in theory is expensive. Thus, we adopt curve-fitting method. Instead of simply solving the relationship between dose and development depth, a fitting method combining convolution model is proposed to improve the fitting accuracy.

We must design such an experiment to obtain the photoresist function in one step. That is, lithography equipment performs the forward process (Fig. 2(c)). By measuring the developed morphology of structure, we can obtain function ${M(x,\;y)}$, and then we can establish the mapping relationship (p function) between the aerial energy distribution $\textrm{E}(x,\;y)$ and function ${M(x,\;y)}$. Subsequently, in the future optimized experiment, the obtained mapping curve can represent the properties of this thick-film photoresist under the fixed lithography process conditions. Aerial energy $\textrm{E}(x,\;y)$ is obtained by convolution of dose distribution ${D(x,\;y)}$ and light source $\textrm{G}(x,\;y)$, aerial energy and dose distribution are in a proportional relationship in theory because convolution is a linear operator. But light source ${G}(x,\;y)$ plays a role of low-pass filtering in this system, and the high-frequency signal of developed pattern will be lost. Therefore, although the middle area of morphology can be optimized to some extent by a linear compensation method, the size of light source must be considered at the edge of the pattern which has more high-frequency signals. In this case, linear compensation can't give an effective solution.

In the LDWL equipment, the input gray-scale image was discretized into a pixel matrix to represent the normalized target structure, and the inverted value of this pixel matrix is operated as the normalized exposure dose matrix. Therefore, by changing the dose matrix and combining with the photoresist characteristic curve fitted by the pre-experiment, we can control the depth and morphology of the developed structure in the photoresist.

2.3 Optimizer

The optimization task is to determine a dose distribution ${D(x,\;y)}$ to change the aerial image distribution $\textrm{E}(x,\;y)$ so that the developed morphology ${M(x,\;y)}$ is closer to the target structure $T{(x,\;y)}$. This optimization task can be expressed by a double integral in Eq. (4). Here, F is defined as the cost function, $\tilde{M}{(x,\;y)}$ represents the estimate of the developed morphology ${M(x,\;y)}$ from Eq. (3), $T{(x,\;y)}$ is the target morphology, and S is our workspace, with width of Nx and height of Ny. Considering the processing of the input grayscale image by the equipment: the value of ${D(x,\;y)}$ is constrained in the range of 0–1, which can be converted into an unconstrained problem using a trigonometric function in Eq. (5).

Combining Eqs. (1)–(5), convolution differential property and convolution theorem which is provided in the Appendix A, derivation of the cost function F in Eq. (4) to in Eq. (5) which transformed convolution into multiplication by Fourier transform (FT) and inverse FT is expressed in Eq. (6). However, the size of light source and workspace are limited, and some high-frequency information in Fourier frequency domain is lost. We can only obtain the approximation of Eq. (7). Therefore, it is important to accurately simulate the light source.

To adapt the optimizer to the pixelated model of LDWL, the above continuous function space is discretized to obtain the discrete matrix space, and the cost function F in Eq. (4) can be expressed in matrix given by Eq. (7), where $\left\| \cdot \right\| _\textrm{2}^\textrm{2}$ is the square of L₂-norm. Matrixes G, D, E, M, T, $\boldsymbol{\omega}$, and F represent ${G(x,\;y)}$, ${D(x,\;y)}$, $\textrm{E}(x,\;y)$, ${M(x,\;y)}$, $T{(x,\;y)}$, $\mathrm{\omega} (x,\;y)$, and F, respectively, after discretization, and variable k represents the number of iterations.

Correspondingly, Eq. (6) in the discrete space is expressed in Eq. (8), and $ \odot $ is the element-by-element multiplication operator. The continuous Fourier transform is replaced by fast Fourier transform (FFT) to accelerate the iteration.

To solve this optimization task reliably, we need to adopt an appropriate optimizer based on derivative. Recently, a constrained gradient search method (L-BFGS) has been applied to the large-scale 2D dose optimization problem to reduce the memory occupation and has shown excellent performance [21]. But different from binary pattern optimization task, morphology of 3D micro/nano structures changes continuously at the edge and is very sensitive to the gradient, which requires us to be more careful when optimizing parameters to avoid the cost function falling into a local minimum.

Here, we are going to try a modified gradient descent method—(Adaptive Moment Estimation) Adam [24]. The optimization results will be compared with L-BFGS, which shows that the optimization error is smaller. We think the reason is that Adam method uses exponential moving average method to modify the iteration step size and direction, so the step size near the local minimum is still enough to escape from the “trap” in such a complicated nonconvex optimization problem.

Adam method has the form in Eq. (9), where ${{\boldsymbol V}_k}$ and ${{\boldsymbol S}_k}$ represent exponential moving average direction vector and exponential moving average step size factor in the k step; and ${\beta _1}$., ${\beta _2}$ represent iterative step size and attenuation coefficient, respectively; $\varepsilon $ considers a small value; and $\mathrm{\alpha}$. is expanded to 1000 times of the common value of 0.003 to adapt to the range of $\mathrm{\omega}$ values. Adam's idea is to use the first and second moments of the gradient combined with moving average method to estimate the descending direction and step size. So the direction and magnitude of the gradient depend not only on the local derivative, but also on the previous path. When the local value of the derivative is 0, there is still a certain “momentum”. The advantage of this is that Adam has a good escape performance for saddle points and some local minima.

Combined with the Adam gradient descent method, numerical iterative optimization can be expressed in the algorithm form in Table 1. The initial dose distribution matrix D₀ is assigned with a matrix $\textrm{0}\textrm{.99}\cdot ({1 - {\boldsymbol T}({x,y} )} )$ attempting to reduce the initialization error, while 1 in D is eliminated to ensure that all values can be iteratively updated in $\mathrm{\omega}$.

Table 1. Optimization algorithm

View Table

3. 3D lithography optimization example

3.1 Target structure

To verify the 3D LDWL model and thick-film photoresist characteristics fitting method of the 3D OPC method and to test the performance of the optimizer in Table 1, we designed and fabricated a group of microlens array (MLA) with a hyperboloid, which has a better aberration elimination effect than that of the spherical lens [26]. The expression of the microlens morphology is presentedn Eq. (9). Conical coefficient is $\textrm{ - 2}\textrm{.25}$, while the radius of the curvature is approximately 150 μm. A single microlens is 80 μm in diameter and 5.25 μm in height.

3.2 Numerical result

The characteristics of the thick-film photoresist are required for the optimization of dose distribution, which often must be obtained through a pre-experiment. The fitting method is provided in Section 2. However, even with the same dose distribution, the change in the experimental conditions will lead to the change of morphology, thereby affecting the photoresist function p, such as photoresist type, developing time, and developer concentration. Therefore, the following experimental parameters of the optimized dose distribution test experiment must be the same as those of the previous experiment. A set of lens grayscale images before 3D OPC were input into LDWL. After exposure and development, a group of lens structures were obtained on the photoresist. Before optimization, these lenses were used for fitting the photoresist function curve (Fig. 3(a)). The normalized dose distribution D corresponding to the gray scale image and normalized light source distribution G measured at the focal plane of the light source were convolved to obtain energy distribution E. The aerial image energy E and development depth at the same position are marked in Fig. 3(a) and fitted well. Positive photoresist was used in the study, which means that higher the energy, the more photoresist is etched.

 

Fig. 3. a) Photoresist function curve (purple line) and measurement data (red dots). The measurement data are obtained from the pre-experiment. X-axis is the convolution result of normalized dose distribution D and normalized light source distribution G, while Y-axis is the etching depth corresponding to the aerial energy E. b) The error descending curve of cost function F is expressed in Eq. (7) (optimized by the Adam gradient descent method). c) Comparison of minimum error between L-BFGS and Adam method. Microlenses with different diameters are used for optimization, and Adam method always has a smaller minimum error in optimization tasks.

Download Full Size | PDF

When the photoresist function p was obtained, the optimization process of dose distribution D made the output of the nonlinear response system composed of LDWL and photoresist closer to the target value T. The measurement shows that the maximum single exposure depth was approximately 7.5 μm, and the target lens was placed at 1 μm above the maximum single exposure depth, with the structure height of 5.25 μm and the diameter of 80 μm. T was normalized when it was input in the initial matrix D₀. MATLAB R2019a was used for numerical simulation and optimization. RTX 1050Ti GPU and Intel Core i7-10700 CPU were used. Every 50 times, a set of data was saved, and for a total of 750 times, optimization consumed 124 s when the grid size was 100 nm/pixel. The error descending curve of cost function F is shown in Fig. 3(b). In addition, microlenses with different sizes ranging from 30 μm to 300 μm in diameter are used for numerical optimization test between Adam and L-BFGS to show the advantages of Adam method in finding the global minimum. The results are given in Fig. 3(c).

To understand the optimization process of Adam, the dose distribution D, aerial energy E and simulated morphology M of 80 μm microlens in different optimization stages are given in Fig. 4. When iteration number ${k\; = 1}$, the input initial dose distribution D₀ which was assigned with a matrix $\textrm{0}\textrm{.99}\cdot ({1 - {\boldsymbol T}({x,y} )} )$ was similar to the inverted T curve. When ${k\; = 300}$, the error was reduced by more than one order of magnitude, but an “oscillation” curve appeared in the top area of the lens morphology in Fig. 4. After 750 times of iteration, the simulation result was close to the target shape, while the error reduced to 2 × 10³(1997). As we can observe from the descending curve, cost function F dropped rapidly before ${k\; = 200}$ times, and then a slight “overshoot” occurred from ${k\; = 200}$ to ${k\; = 400}$ times. After 400 times, the curve became horizontal. We speculate that “overshoot” in the descending curve was related to the appearance “oscillation” of morphology M. One possible explanation is that the morphology error values of different regions o. the microlens in F were different. When the step size was suitable for the surrounding area, it was too large for the top area, resulting in “oscillation”. But it will always disappear with the increase of iteration number k. Although there is a short process of error increase in the optimization by Adam method, it gives Adam more opportunities to adjust the direction to find the global minimum. The comparison in Fig. 3(c) shows that Adam always has a smaller minimum error in the morphology correction of micro/nano structures than L-BFGS method.

 

Fig. 4. Iterative process of the optimization; the number of iterations increases from left to right. This shows the side views of dose D, aerial energy E, and sulated morphology M at $k = 1,\; 300\; and\; 750$ times (given in different colors). The simulated morphology results at$\; k\; = \; 1\; $are close to the experimental results shown in Fig. 5(e). When $k = 300$, the morphology is closer to the target but “overshoot.” When $k = 750$, the error between the optimized morphology and the target is reduced to a small level (coincident in the figure). In this process, the loss descending curve of the cost function in Eq. (7) decreases rapidly, which is shown in Fig. 3(b). All the unit of the X-axis in the figure is in 100 nm/pixel.

Download Full Size | PDF

3.3 Numerical result discussion

At present, we have discussed the stability and global optimization effect of Adam optimizer. The stability and efficiency of the optimizer are acceptable, and the error from morphology can converge to a relatively small value in optimization. But we haven't verified the reliability of the 3D LDWL model and the optimization result of the optimizer. These verifications must be carried out by lithography experiments.

4. Experiment

In the pre-experiment and verification experiment, the LDWL equipment was Picomaster-100 (4PICO Litho Co., Ltd). The spot size (FWHM) was 850 nm with 200 nm step. The maximum exposure dose was 1000 mj/cm², while the exposure threshold was 50 mj/cm². We used a 10-inch photoresist plate made by Hunan Puzhao Co., Ltd. Photoresist AZ4562 was evenly applied to the glass plate. Finally, structures were developed in 20% KOH at 22 ℃ for 140 s.

4.1 Morphological error

The measurement results from optical profiler (ZYGO Co., Ltd, Nexview NX2) of MLA before and after 3D OPC in the subsequent experiment are shown in Fig. 5. The observation results of the microlens before and after 3D OPC under microscope are shown in Figs. 5(a) and 5(c), respectively. The 3D physical appearance captured by optical profiler is shown in Figs. 5(b) and 5(d). The dose range of MLA before 3D OPC was adjusted to minimize the error of morphology and target, and the purpose was to highlight the errors caused by the morphology before optimization. Significant differences can be observed in physical appearance and 3D view. Figures 5(e) and 5(f) show the respective measurement results. As shown in Fig. 5 g, compared with the 0.955 μm (18.2% in height) average absolute error of the before 3D OPC structure, the 0.274 μm (5.22% in height) average absolute error of the optimized structure was smaller, which is closer to the errors possibly caused by process and measurement.

 

Fig. 5. Physical appearance of the lens before a) and after c) 3D OPC under a 20x microscope. 3D view of MLA before b) and after d) 3D OPC under observation of optical profiler, respectively. e) Morphology M before (in red) and after (in blue) 3D OPC compared with the target morphology (in black). f) The absolute height error between the lens and target at different positions before (in red) and after (in blue) 3D OPC. The horizontal lines show the average error values before (in red) and after (in blue) 3D OPC.

Download Full Size | PDF

4.2 Image performance of MLA

For industrial applications, comparing the imaging effect after transfer printing is necessary, because the ultimate goal is always to improve the performance of the micro/nano structures. Polydimethylsiloxane (PDMS) is widely used for batch replication of micro-optical devices because of its good light transmission performance and transfer precision of micro/nano structures [27–29]. Here, we transferred the structure to PDMS and compared the imaging effect of MLA before and after 3D OPC. The imaging object was the 1951 U.S. Air Force resolution test card, which had eight groups of patterns from 0 to 7, and each group had six groups of mutually perpendicular stripe patterns. Figure 6 shows the image of portion of the test card under MLA before and after 3D OPC. Compared with the structure before 3D OPC, the field of view (FOV) was expanded from 11.7% to 57.0% in the area (from 34.2% to 75.5% in the width and height). The increase of the FOV intuitively shows the improvement of imaging performance due to the optimized morphology. Owing to space limitations, only the improvement of the FOV of microlens after 3D OPC is provided here. For more detailed measurement of microlens imaging, please refer to our related work [30].

 

Fig. 6. Resolution card (partial. imaged by MLA before a) and after b) 3D OPC and resolution card (partial, as a control). The fields of view (FOV, in white block) account for 11.7% and 57.0% of the total area (in red block) before and after 3D OPC, respectively (34.2% and 75.5% in width and height, respectively). Correspondingly, the imaging is clearer.

Download Full Size | PDF

4.3 Experiment discussion

In this section, we apply the optimized dose distribution to a lithography experiment. The morphology error and final imaging effect in the experiment are greatly improved compared with those before 3D LDWL. In the experiment of the microlens array, the manufacturing error of the LDWL device was limited to 5.22%, and the lens morphology was close to our expectation. However, manufacturing errors still exist. Owing to the combination of the incoherent 3D LDWL model and thick-film photoresist, errors may emerge from two aspects: the LDWL model is not accurate enough, such as errors caused by the simplification of the light source model, and the beam is diffracted to a wide area, and diameter of the light source model is limited, which has been smoothed by Gaussian. The process error, for example, the photoresist surface, cannot be completely horizontal, and measurement error may be imported by the fitting curve. In addition, the structure with strict requirements for the overall height and placement position can lead to more requirements for process repeatability.

Based on the above analysis, in order to further improve the optimization accuracy, it is necessary to consider more process parameters. The 3D OPC model proposed by the study has the expansibility of process parameters. The key point is that the process parameters are related to the photoresist characteristic curve, and different process parameters have different photoresist curves. For example, if the development time is changed, the development depth will also change, and the photoresist characteristic curve will be changed. Therefore, the optimal process parameters can be determined only by comparing the optimized errors of different process parameters. A general flow is shown in Fig. 7.

 

Fig. 7. Extended 3D OPC method considering process parameters. Different process parameters correspond to different photoresist characteristic curves. 3D OPC method optimizes the dose distribution of different process parameters and compares them to get the optimal dose distribution and corresponding process parameters, such as developing time and solvent ratio of developer.

Download Full Size | PDF

In our additional experiments, 3D OPC method was applied to other structures, such as slope, stairs, and sidewall to further prove the reliability. Please refer to Appendix B, Fig. 8, Fig. 9, and Fig. 10.

5. Conclusion

LDWL has great advantages in the fabrication of high-degree-of-freedom morphology elements; however, it is has limitations owing to the influence of many nonlinear factors. In this study, combining the point-by-point exposure characteristic of LDWL with Fourier analysis, a 3D OPC method was proposed to correct the distortion of the LDWL fabrication structures. The model was established in the second quarter, applied in the third quarter, and verified through experiments in the fourth quarter. The average lens morphology error was controlled within 5.22%, which is often acceptable. The error comparison results are shown in Fig. 5(f)) between the pre-and post-optimization morphology and the target morphology. The improvement of the FOV of the optimized lens show the effectiveness of the method. In our additional experiments, 3D OPC method was applied to other structures, such as slope, stairs, and sidewall. The results indicate that 3D OPC can be widely used in micro/nanostructures with various morphologies and has a wide application prospect in micro/nano optics, silicon optical chips, microfluidic, superhydrophobic structure, and MEMS. Using 3D OPC to correct the morphology of micro/nanostructures can better realize their function. The 3D OPC method combining the photoresist function to optimize the exposure dose of a thick-film photoresist provides an acceptable approach for the LDWL system to improve the manufacturing accuracy and reduced the production and experiment cost. In the future work, we will consider more process parameters in the 3D OPC model to further improve the optimization accuracy. Simultaneously, we will also study the resolution limit of smaller structures with a certain light source and more possibilities of 3D OPC applications.

Appendix A: derivation of gradient ${{\partial {\boldsymbol F}} / {\partial {\boldsymbol \omega }}}$

Differential properties of convolution:

(11)$$\begin{array}{l} \frac{{\partial (D \otimes G)}}{{\partial \omega }}\\ \textrm{ } = \frac{\partial }{{\partial \omega }}(\int_0^{Nu} {\int_0^{Nv} {(D(u,v)) \cdot G(x - u,y - v)dudv} } )\\ by\textrm{ }Leibniz\textrm{ }formula:\\ \textrm{ } = \int_0^{Nu} {\int_0^{Nv} {\frac{\partial }{{\partial \omega }}((D(u,v)) \cdot G(x - u,y - v))dudv} } \\ \textrm{ } = \int_0^{Nu} {\int_0^{Nv} {\frac{\partial }{{\partial \omega }}(D(u,v)) \cdot G(x - u,y - v)dudv} } \\ \textrm{ } = \frac{{\partial D}}{{\partial \omega }} \otimes G\\ \textrm{ } = \frac{{\partial E}}{{\partial \omega }} \end{array}$$

(12)$$\begin{array}{l} \int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {\sin \omega (u,v) \cdot G(x - u,y - v)dudv} } \\ \textrm{ } = F{T^{ - 1}}(FT(\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {\sin\omega (u,v) \cdot G(x - u,y - v)dudv} } ))\\ \textrm{ } = \textrm{F}{\textrm{T}^{ - 1}}(\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {\sin\omega (u,v) \cdot G(x - u,y - v)dudv} } } } {e^{ik(x + y)}}dpdq)\\ let\textrm{ }p = x - u,q = y - v:\\ \textrm{ } = \textrm{F}{\textrm{T}^{ - 1}}(\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {\sin\omega (u,v) \cdot G(p,q)dudv} } } } {e^{ik(p + u + q + v)}}dpdq)\\ \textrm{ } = \textrm{F}{\textrm{T}^{ - 1}}(\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {\sin\omega (u,v){e^{ik(u + v)}}dudv \cdot \int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {G(p,q)} } } } {e^{ik(p + q)}}dpdq)\\ \textrm{ } = \textrm{F}{\textrm{T}^{ - 1}}(FT(\sin\omega (u,v)) \cdot FT(G(p,q))) \end{array}$$

Equation (11) is accelerated by convolution theorem, Fourier transform (FT) is used to derive the formula, and fast Fourier transform (FFT) replaced FT in a computer to accelerate.

Combined Eqs. (1)–(5), (10), and (11). Equation (6) is derived as follows:

(13)$$\begin{array}{l} \frac{{\partial F}}{{\partial \omega }} = \textrm{2}\int_0^{Nx} {\int_0^{Ny} {(\tilde{M} - T) \cdot \frac{{\partial \tilde{M}}}{{\partial \omega }}dxdy} } \\ \textrm{ } = \textrm{2}\int_0^{Nx} {\int_0^{Ny} {(\tilde{M} - T) \cdot \frac{{\partial p(E)}}{{\partial \omega }}dxdy} } \\ \textrm{ } = \textrm{2} \cdot \dot{p}(E) \cdot \int_0^{Nx} {\int_0^{Ny} {(\tilde{M} - T) \cdot \frac{{\partial E}}{{\partial \omega }}dxdy} } \\ \textrm{ } = \textrm{2} \cdot \dot{p}(E) \cdot \int_0^{Nx} {\int_0^{Ny} {(\tilde{M} - T) \cdot (\frac{{\partial D}}{{\partial \omega }} \otimes G)dxdy} } \\ \textrm{ } = \textrm{2} \cdot \dot{p}(E) \cdot \int_0^{Nx} {\int_0^{Ny} {(\tilde{M} - T) \cdot ( - \frac{1}{2}\sin \omega \otimes G)dxdy} } \\ \textrm{ } = \dot{p}(E) \cdot \int_0^{Nx} {\int_0^{Ny} {(T - \tilde{M}) \cdot (\sin \omega \otimes G)dxdy} } \\ \textrm{ } = \dot{p}(E) \cdot \int_0^{Nx} {\int_0^{Ny} {(T - \tilde{M}) \cdot \textrm{F}{\textrm{T}^{ - 1}}(\textrm{FT}(\sin \omega ) \cdot \textrm{FT}(G))dxdy} } \end{array}$$

Appendix B: more experiments

To ensure the reliability of the 3D OPC method. We carried out experimental verification on more structures with different morphology. For the slope, the average absolute error is improved from 0.821μm (20.5%) to 0.191μm (4.78%), while the average absolute error of stairs is optimized from 0.756μm (18.9%) to 0.214μm (5.34%).The FWHM and sidewall angle are also improved (2.70μm (67.5% in total 4μm) to 3.03μm (75.8% in total 4μm), 78° to 80°).

 

Fig. 8. Slope structure is tested. a)c) and b)d) are the microscope observation results and the 3D morphology before (a),b)) and after (c),d)) 3D OPC. e)f) showed the target morphology (in black) and the measurement results before and after optimization. The results show that the average absolute error of slope morphology is optimized from 0.821µm (20.5%) to 0.191µm (4.78%).

Download Full Size | PDF

 

Fig. 9. Stairs structure is tested. The results show that the average absolute error of stairs morphology is optimized from 0.756µm (18.9%) to 0.214µm (5.34%).

Download Full Size | PDF

 

Fig. 10. Sidewall structure is tested. The FWHM and sidewall angle are both improved (2.70µm (67.5% in total 4µm) to 3.03µm (75.8% in total 4µm), 78° to 80°).

Download Full Size | PDF

Funding

National Natural Science Foundation of China (U20A6004).

Acknowledgments

This work was supported by the National Natural Science Foundation of China (NO. U20A6004). Thanks for the Simax Shanghai Co., Ltd technical support.

Disclosures

There are no financial conflicts of interest to disclose.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. N. Tsutsumi, J. Hirota, K. Kinashi, and W. Sakai, “Direct laser writing for micro-optical devices using a negative photoresist,” Opt. Express 25(25), 31539 (2017). [CrossRef]  

2. Z. Tian, W. Yao, J. Xu, Y. Yu, Q. Chen, and H. Sun, “Focal varying microlens array,” Opt. Lett. 40(18), 4222 (2015). [CrossRef]  

3. A. Greiner, B. Richter, and M. Bastmeyer, “Micro- and nano-engineered 3D scaffolds for cell culture studies,” Macromol. Biosci. 12(10), 1301–1314 (2012). [CrossRef]  

4. P. Brenner, O. Bar-On, T. Siegle, T. Leonhard, R. Gvishi, C. Eschenbaum, H. Kalt, J. Scheuer, and U. Lemmer, “3D whispering-gallery-mode microlasers by direct laser writing and subsequent soft nanoimprint lithography,” Appl. Opt. 56(13), 3703–3708 (2017). [CrossRef]  

5. L. Gao, X. Wang, W. Dai, R. Bagheri, and G. Zou, “Laser direct writing assisted fabrication of skin compatible metal electrodes,” Adv. Mater Technol. 5(5), 2000012 (2020). [CrossRef]  

6. Y. Cheng, Y. Liao, F. He, L. Qiao, and Z. Xu, “Fabrication of microfluidic systems with arbitrary 3D geometries inside fused silica using femtosecond laser direct writing,” in International Photonics and Optoelectronics Meetings (Optica Publishing Group, Wuhan, 2012), pp. h4C–h5C.

7. S. Luan, P. Xu, Y. Zhang, L. Xue, Y. Song, and C. Gui, “Flexible Superhydrophobic Microlens Arrays for Humid Outdoor Environment Applications,” ACS Appl. Mater. Inter. 14(47), 53433–53441 (2022). [CrossRef]  

8. A. J. Fleming, A. G. Wills, and B. S. Routley, “Exposure optimization in scanning laser lithography,” IEEE Potentials 35(4), 33–39 (2016). [CrossRef]  

9. Y. Hirai, K. Sugano, T. Tsuchiya, and O. Tabata, “A three-dimensional microstructuring technique exploiting the positive photoresist property,” J. Micromech Microeng. 20(6), 065005 (2010). [CrossRef]  

10. W. T. Plummer, J. G. Baker, and J. Van Tassell, “Photographic optical systems with nonrotational aspheric surfaces,” Appl. Opt. 38(16), 3572–3592 (1999). [CrossRef]  

11. L. Yang, F. Shen, Z. Ding, X. Tao, Z. Zheng, F. Wu, Y. Li, and R. Wu, “Freeform optical design of beam shaping systems with variable illumination properties,” Opt. Express 29(20), 31993–32005 (2021). [CrossRef]  

12. T. Xionggui, G. Yongkang, and D. Jingleiet, “Analysis of diffractive image field in thick film photoresist by using Fouriermodal method,” Acta Optica Sinica, 246–250 (2005).

13. F. C. M. van Kempen, Y. Hirai, F. van Keulen, and O. Tabata, “Automatic process design for 3D thick-film grayscale photolithography,” in Proceedings of Transducers & Eurosensors XXVII (IEEE, 2013), pp. 1625–1628.

14. A. Erdmann, C. L. Henderson, and C. G. Willson, “Impact of exposure induced refractive index changes of photoresists on the photolithographic process,” J. Appl. Phys. 89(12), 8163–8168 (2001). [CrossRef]  

15. Y. Sensu and A. Sekiguchi, “Improved resolution of thick film resist (effect of pre-bake condition),” Proc. SPIE 4979, 561–581 (2003). [CrossRef]  

16. G. E. Flores, W. W. Flack, and E. Tai, “Investigation of the properties of thick photoresist films,” Proc. SPIE 2195, 734–751 (1994). [CrossRef]  

17. J. Dunkel, F. Wippermann, A. Brückner, A. Bräuer, and A. Tünnermann, “Laser lithographic approach to micro-optical freeform elements with extremely large sag heights,” Opt. Express 20(4), 4763 (2012). [CrossRef]  

18. X. Xiao, Study on imaging characteristics of thick-film resist, Chap. 4, (Sichuan University, 2003).

19. A. J. Fleming, A. Wills, O. T. Ghalehbeygi, B. Routley, and B. Ninness, “A nonlinear programming approach to exposure optimization in scanning laser lithography,” in Proceedings of 2016 American Control Conference (ACC) (American Automatic Control Council (AACC), 2016), pp. 5811–5816.

20. A. Poonawala and P. Milanfar, “Mask design for optical microlithography—An inverse imaging problem,” IEEE Trans. on Image Process. 16(3), 774–788 (2007). [CrossRef]  

21. C. Jidling, A. J. Fleming, A. G. Wills, and T. B. Schön, “Memory efficient constrained optimization of scanning-beam lithography,” Opt. Express 30(12), 20564–20579 (2022). [CrossRef]  

22. J. Li, S. Liu, and E. Y. Lam, “Efficient source and mask optimization with augmented Lagrangian methods in optical lithography,” Opt. Express 21(7), 8076 (2013). [CrossRef]  

23. X. Sun, S. Yin, H. Jiang, W. Zhang, M. Gao, J. Du, and C. Du, “U-Net convolutional neural network-based modification method for precise fabrication of three-dimensional microstructures using laser direct writing lithography,” Opt. Express 29(4), 6236 (2021). [CrossRef]  

24. D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” CoRR abs/1412.6980 (2015).

25. T. Xionggui, G. Fuhua, and G. Fenget, “Study on distribution of PAC concentration in thick film resist,” Microfabricaion Technology, 48–53 (2005).

26. E. Heynacher, “Aspheric optics. How they are made and why they are needed,” Phys. Technol. 10(3), 124–131 (1979). [CrossRef]  

27. Y. Shu, T. Chen, W. Zhou, Z. Zhou, and A. Y. Yi, “Fabrication of large-scale infrared diffractive lens arrays on chalcogenide glass by means of step-and-repeat hot imprinting and non-isothermal glass molding,” Int. J. Adv. Manuf. Technol. 116(9-10), 3075–3085 (2021). [CrossRef]  

28. P. Chen, R. Zhang, and L. Chen, “Using micromachined molds, partial-curing PDMS bonding technique, and multiple casting to create hybrid microfluidic chip for microlens array,” Micromachines 10(9), 572 (2019). [CrossRef]  

29. Q. Li, M. G. Ji, and J. Kim, “Grayscale nanopixel printing at sub-10-nanometer vertical resolution via light-controlled nanocapillarity,” ACS Nano 14(5), 6058–6066 (2020). [CrossRef]  

30. S. Luan, F. Peng, G. Zheng, C. Gui, Y. Song, and S. Liu, “High-speed, large-area and high-precision fabrication of aspheric micro-lens array based on 12-bit direct laser writing lithography,” Light: Adv. Manufac. 3, 1 (2022). [CrossRef]