Topology Optimization Design of an Active Deformable Mirror Based on Discrete Orthogonal Zernike Polynomials

Abstract

In order to design an active deformation mirror for projection objective aberration imaging quality control, a topology optimization design method of active deformation mirrors based on discrete orthogonal Zernike polynomials is proposed in this paper. Firstly, in order to solve the problem that continuous Zernike polynomials do not have orthogonality on the discrete coordinates inside the unit circle, which causes the instability of topology optimization results, discrete orthogonal Zernike polynomials are used to characterize the active deformation mirror wave aberrations. Then, the optical and structural deformations are combined to establish an optical-mechanical coupling topology optimization model with the help of the variable density method to derive the sensitivity of the mathematical model. Finally, a wave aberration corrected deformation mirror in an optical machine system is used as an arithmetic example for topology optimization, and the results show that the absolute value of the Zernike coefficient Z₄ after optimization is improved by nearly one order of magnitude compared with the value before optimization, and the vibration characteristics of the optimized structure meet the design requirements. The optimization effect is significant, which improves the optical performance of the deformed mirror and provides a new scheme for the design of the deformed mirror structure which has a certain practical value for engineering.

1. Introduction

Lithography is the core technology of very large-scale IC manufacturing. With the development of integrated circuit transistors of 22 nm, 10 nm, 7 nm or even smaller sizes, the traditional immersion lithography has gradually failed to meet the production requirements and the cost of other lithography technologies such as double lithography is rising rapidly which makes the professionals in the industry once again expect extreme ultraviolet lithography (EUVL) to support the continuation of Moore’s law in semiconductor manufacturing, and EUVL is recognized as the technology with the most potential to achieve integrated circuit manufacturing technology to enter nodes below 1 nm [1,2].

The resolution of the projection objective lens, the core component of the EUVL mask aligner, is increasing higher and higher, which constantly challenges the limits of the existing technologies making it is urgent to develop new technologies. The resolution (RES) of the projection objective is related to the working wavelength $\lambda$, numerical aperture NA, and process factor k₁, which is expressed as

$$RES=k_1\frac{\lambda }{NA}$$

(1)

Equation (1) guides the development direction of lithographic projection objectives of shorter working wavelengths, larger numerical apertures, and smaller process factors. The development route of decreasing wavelengths and increasing numerical apertures has brought great challenges to the design and manufacture of the projection objective, the core component of a mask aligner [3].

Limited by the diffraction effect, the existence of wave image aberration can always cause the redistribution of light from the center of the image to the diffraction ring. It usually causes edge rounding, distortion, line length shortening, line width deviation of sparse and dense lines, and line width differences between transparent and opaque masks in spatial images, which seriously affects the resolution of the EUVL projection system.

In order to control the objective aberration and make the imaging quality approach the diffraction limit, the active deformation mirror is used as the wavefront correction device to compensate the wavefront aberration in the projection objective lens, which has become a problem that must be solved by researchers of lithography objectives. Using structural optimization design methods to design active deformation mirrors including size optimization, shape optimization, and topology optimization, has become an important current research topic [4]. Size optimization and shape optimization methods are easy to operate and are widely used, however they can only change the size, shape, and morphology of the active deformation mirror, and cannot optimize the topological configuration of the active deformable mirror structure, so the optimization benefit is low. Topology optimization, by designing the shape, number, and distribution position of the internal holes of the structure, obtains the optimal topological configuration of the structure, while optimizing the shape of the external boundary of the structure. Its biggest advantage is that under the premise that the specific structure is unknown, according to the known boundary conditions and load conditions, the optimal topology can be obtained after optimization, which also provides designers with new design schemes and optimal material allocation schemes.

Topology optimization is one of the most advanced techniques because it has many variables and can remove unnecessary units. Common methods include the homogenization method, variable density method, evolutionary algorithm ESO, and level set method. The homogenization method realizes the addition and deletion of microstructures through the growth and decline of microstructure size, and finally obtains the macroscopic topological configuration. It is suitable for material configuration optimization design, but the sensitivity of solving macroscopic equivalent elastic tensor and structural response function is complicated [5,6,7]. The variable density method is inspired by the homogenization method. It establishes the relationship between density and elastic modulus, takes the material density as the topological optimization design variable, and obtains the optimal distribution of structural materials through the presence or absence of materials. Compared with the homogenization method, the variable density method can directly obtain the relationship between the unit density and the elastic modulus without solving it again, which greatly reduces the variables of optimization design, simplifies the optimization solution process, and improves the calculation efficiency [8,9,10,11,12]. The evolutionary ESO method realizes structural topology optimization by imitating the biological evolution process and gradually deleting invalid or inefficient materials. The evolutionary ESO method avoids using mathematical programming to solve large-scale design variables, which greatly improves the efficiency of the optimization solution. Compared with the variable density method, the ESO method lacks strict mathematical theory proof, and its universality, robustness, and numerical stability are inferior to those of the variable density method. Usually, reasonable parameters need to be set in order to obtain satisfactory results [13,14]. The level set method solves the Hamilton–Jacobi partial differential equation to obtain the level set function of the extended domain, and then obtains the topological configuration of the structure. The advantages of the level set method are that in the process of structural topology optimization, the topological changes of geometric boundaries of materials can be handled naturally, the checkerboard effect can be overcome, the configuration of topology optimization results is clear, there are no gray units, and the finite element mesh and optimization solution are separated. However, each iteration needs to solve the complex Hamilton–Jacobi differential equation by difference, which leads to low-computational efficiency and complicated optimization [15,16,17].

Scholars at home and abroad have done substantial research on the application of the topology optimization method in lightweight optical systems [18,19,20,21,22]. Li Chengliang et al. [15] took the surface RMS value as the optimization objective and adopted the level set method to carry out the topology design of the large-aperture mirror. Park et al. use Zernike polynomials to describe the mirror surface shape and takes the mean square difference between the actual and ideal deformation of the mirror surface as the objective function to develop a topology optimization model for the main reflector support structure of a multispectral camera system [23]. The topology optimization results of Park et al. vary with the volume constraints and the topological configuration of the reflector support structure, and the topology optimization result is an unstable problem. Preliminary analysis of the problem is that it may be caused by the Zernike fitting equation constructed directly with finite element discrete nodes, ignoring the fact that the Zernike polynomials do not have the nature of orthogonality in the discrete coordinates inside the unit circle, resulting in a serious pathology in the condition number of the equation coefficient matrix. Small changes in optimization variables will lead to large changes in the Zernike polynomial fitting coefficient, which will lead to large changes in the optimal value of objective function. Finally, the results of the topology optimization show that the topology optimization is unstable. Sahu et al. used the Zernike polynomial coefficient as optimization constraints and HyperWorks to realize topology optimization of a 1.2 m space mirror support structure [24]. Sahu et al. consider the influence of structural design on optical aberration from the perspective of constraint. There is a big error of constraint conditions in the optimization process, which is still unable to meet the requirements of imaging quality for wave aberration. Koppen et al. carried out research on the topology optimization method of two-dimensional mirror support structures [25]. In their study, the mean square deviation of actual deformation and ideal deformation was taken as the objective function, and the relationship between the optimization objective function and the surface shape error was expounded. However, the topology optimization problem with wave aberration as the optimization objective was still not solved in the objective function, and this method only optimizes the two-dimensional plane model. In addition, these studies only used the topology optimization results as a reference for the layout and lightweight design of stiffeners.

In order to design active deformation mirrors suitable for aberration imaging quality control of projection objective, this paper proposes an active deformation mirror topology optimization design method based on optical force coupling. Firstly, to address the problem that the Zernike polynomials are not orthogonal in the discrete coordinates inside the unit circle and the small changes of the optimization variables will lead to large changes of the wave aberration and thus cause the topological optimization results to be unstable, the discrete orthogonal Zernike polynomials are used to characterize the wave aberration of the active deformation mirror. Then, the optical and structural deformation are combined, and the optical and mechanical coupling topology optimization model is established using the variable density method, and the sensitivity of the mathematical model is deduced. Finally, the wave aberration correction deformable mirror in an opto-mechanical system is taken as an example, and the optimal topology configuration of the deformable mirror satisfying the convergence condition is obtained based on the variable density algorithm under the weight constraint condition, which provides conceptual design and reference for the opto-mechanical structure design.

2. Mathematical Model for Topology Optimization of Wave Aberration Structure

Optical wave aberration, which is commonly used to assess the quality and the shape of the wave aberration, is rather complicated and in general it is difficult to use its simple function accurately, in all kinds of approximate method by Zernike polynomials in the circle domain orthogonal to each other, and the low-order model and the corresponding seidel aberration in optical design become important to describe the wave aberration of the description method [26].

The relationship between the Zernike polynomial and the Seidel aberration is as follows.

$$\begin{array}{l}Seide{l}_{Tilt.mag}=\sqrt{Z_2^2+Z_3^2}\\ Seide{l}_{Tilt.angle}=a\tan (\frac{Z_3}{Z_2})\\ Seide{l}_{Focus.mag}=2Z_4-6Z_9\\ Seide{l}_{Asti.mag}=2\sqrt{Z_5^2+Z_6^2}\\ Seide{l}_{Asti.angle}=\frac{1}{2}a\tan (\frac{Z_6}{Z_5})\\ Seide{l}_{Coma.mag}=3\sqrt{Z_7^2+Z_8^2}\\ Seide{l}_{Coma.angle}=a\tan (\frac{Z_8}{Z_7})\\ Seide{l}_{Sphe.mag}=6Z_9\end{array}$$

(2)

where, mag represents size, and tilt, focus, Asti, coma, and sphe respectively represent tilt, defocus, primary astigmatism, primary coma, and primary spherical aberration in Seidel aberration. Z₂ and Z₉ are Zernike polynomial coefficients.

The optical element in its support structure, under the action of load and the external environment causes the deformation of the optical element. The optical element deformation control of the imaging quality is described by wave aberration. The deformation of the structure topology optimization technology to structure topology optimization model was set up, to describe wave aberration enable Zernike polynomial coefficients of the objective function to establish the topology optimization model. By linking the image quality control of optical system with structural topology optimization, a mathematical model of wave aberration structure topology optimization is formed. The mathematical model based on the finite element discrete form is as follows:

$$\min f(U,\tilde{\rho })=L_1^T{Z}_{coeff}(U,\tilde{\rho })$$

(3)

$$\mathrm{s}.\mathrm{t}.KU=F$$

(4)

$$Z_{fund}{Z}_{coeff}=L_{2}U$$

(5)

$$V_{\min }\le V=\frac{{\int }_{\mathrm{\Omega }}{\tilde{\rho }}_{e}d\mathrm{\Omega }}{V^{*}}\le V_{\max }$$

(6)

$$0\le \tilde{\rho }\le 1$$

(7)

The specific form of Equation (5) is as follows:

$$\begin{array}{l}{Z}_{fund1}(r_1,{\theta }_0)Z_{coeff1}+Z_{fund2}(r_1,{\theta }_0)Z_{coeff2}+\cdots +Z_{fundN}(r_1,{\theta }_0)Z_{coeff\mathrm{N}}=U(r_1,{\theta }_0)\\ Z_{fund1}(r_i,{\theta }_j)Z_{coeff1}+Z_{fund2}(r_i,{\theta }_j)Z_{coeff2}+\cdots +Z_{fundN}(r_i,{\theta }_j)Z_{coeff\mathrm{N}}=U(r_i,{\theta }_j)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\cdots \hspace{0.17em}\hspace{0.17em}\hspace{1em}\hspace{1em}\\ Z_{fund1}(r_I,{\theta }_J)Z_{coeff1}+Z_{fund2}(r_I,{\theta }_J)Z_{coeff2}+\cdots +Z_{fundN}(r_I,{\theta }_J)Z_{coeff\mathrm{N}}=U(r_I,{\theta }_J)\end{array}$$

(8)

$$Z_{fund}=\left[\begin{array}{l}{Z}_{fund1}(r_1,{\theta }_0)\hspace{0.17em}\hspace{0.17em}{Z}_{fund2}(r_1,{\theta }_0)\cdots Z_{fundN}(r_1,{\theta }_0)\\ Z_{fund1}(r_i,{\theta }_j)\hspace{0.17em}\hspace{0.17em}{Z}_{fund2}(r_i,{\theta }_j)\cdots Z_{fundN}(r_i,{\theta }_j)\\ \hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\cdots \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\cdots \hspace{1em}\hspace{1em}\hspace{0.17em}\cdots \hspace{1em}\hspace{1em}\hspace{1em}\\ Z_{fund1}(r_I,{\theta }_J)\hspace{0.17em}{Z}_{fund2}(r_I,{\theta }_J)\cdots Z_{fundN}(r_I,{\theta }_J)\end{array}\right]$$

(9)

The Zernick polynomial Z(r, θ), which can be resolved as a function of the radial coordinate r and the angular coordinate θ, is as follows.

$$\begin{array}{l}{Z}_{fundk}(r,\theta )=Z_n^m(r,\theta )=\sqrt{n+1}{R}_n^m(r)e^{im\theta }\\ R_n^m(r)={\displaystyle \sum _{s=0}^{(n-m)/2}\frac{{(-1)}^s(n-s)!}{s!\left(\frac{n+m}{2}-s\right)!\left(\frac{n-m}{2}-s\right)!}}{r}^{n-2s}\end{array}$$

(10)

where $Z_{fund}$ is the matrix of the Nth-order Zernike polynomial basis functions and $Z_{coeff}$ is the vector of the Nth-order Zernike polynomial coefficients, which is obtained by Equation (5) using the least-squares method. Equation (4) is the finite element equilibrium equation, K is the overall stiffness matrix of the finite element discrete equation, U is the displacement solution vector of the finite element discrete equation, F is the equivalent nodal load, and Equation (6) is the volume constraint. The vectors L₁ and L₂ are vectors or matrices containing only 0 and 1 to extract the corresponding Zernike polynomial coefficients and displacements, $\tilde{\rho }$ are the optimal design variables and $r$, $\theta$ are polar coordinates.

The interpolation model is based on the improved SIMP method, which has the advantages of simple modelling, convenient solution, high-computational efficiency, universality, and portability. The interpolation model is as follows:

$$E=E_{\min }+{\tilde{\rho }}_e^P(E_0-E_{\min })$$

(11)

where E is the macroscopic modulus of elasticity of the material; E_min is a very small positive value used to ensure the nonsingularity of the overall stiffness matrix; E₀ is the modulus of elasticity of the solid material; p is the penalty factor; and ${\tilde{\rho }}_e$ is the relative physical density of element e. ${\tilde{\rho }}_e=1$ indicates a solid cell and ${\tilde{\rho }}_e=0$ indicates that the cell has no material distribution. Here ${\tilde{\rho }}_e$ is expressed as a function of the nodal design variable ${\tilde{\rho }}_e$.

The Lagrange multiplier method is used to process the constrained optimization model. The Lagrangian function is as follows:

$$L(U,\tilde{\rho })=L_1^T{Z}_{coeff}+{\lambda }_1^T(KU-F)+{\lambda }_2^T(Z_{fund}{Z}_{coeff}-L_{2}U)$$

(12)

The derivatives of the Lagrangian function concerning the optimized design variables are found as follows:

$$\begin{array}{l}\frac{\partial L}{\partial \tilde{\rho }}=L_1^T\frac{\partial Z_{coeff}}{\partial U}\frac{\partial U}{\partial \tilde{\rho }}+{\lambda }_1^T\frac{\partial K}{\partial \tilde{\rho }}U+{\lambda }_1^{T}K\frac{\partial U}{\partial \tilde{\rho }}\hspace{0.17em}+{\lambda }_2^T{Z}_{fund}\frac{\partial Z_{coeff}}{\partial U}\frac{\partial U}{\partial \tilde{\rho }}-{\lambda }_2^T{L}_2\frac{\partial U}{\partial \tilde{\rho }}\\ \hspace{1em}\hspace{0.17em}={\lambda }_1^T\frac{\partial K}{\partial \tilde{\rho }}U+(L_1^T+{\lambda }_2^T{Z}_{fund})\frac{\partial Z_{coeff}}{\partial U}\frac{\partial U}{\partial \tilde{\rho }}+({\lambda }_1^{T}K-{\lambda }_2^T{L}_2)\frac{\partial \widehat{U}}{\partial \tilde{\rho }}\end{array}$$

(13)

Here the adjoint equation is introduced as

$$L_1^T+{\lambda }_2^T{Z}_{fund}=0$$

(14)

$${\lambda }_1^{T}K-{\lambda }_2^T{L}_2=0$$

(15)

The derivative of the Lagrangian function with respect to the optimized design variables is simplified as

$$\frac{\partial L}{\partial \tilde{\rho }}={\lambda }_1^T\frac{\partial K}{\partial \tilde{\rho }}U$$

(16)

In the numerical implementation of the SIMP method, there are two kinds of numerical instability problems: chessboard format and grid dependence. The method of filtering sensitivity and density can effectively avoid numerical instability [27]. In this paper, the weighted filtering method is adopted to filter sensitivity and density. The expression of sensitivity filtering is as follows:

$$\widehat{\frac{\partial f}{\partial {\rho }_j}}=\frac{1}{\max (\gamma ,{\rho }_j){\displaystyle \sum _{i\in N_j}{H}_{ji}}}{\displaystyle \sum _{i\in N_j}{H}_{ji}{\rho }_i}\frac{\partial f}{\partial {\rho }_i}$$

(17)

where f represents the objective function or constraint function, ${\rho }_i$ and ${\rho }_j$ are the densities of node i and node j, respectively, $\gamma$ is a small positive value (generally desirable $\gamma ={10}^{-3}$) to avoid zero values in the denominator, N_j is the set of all neighboring nodes within a radius r_min with node j as the origin, H_ji is the weighting factor, and the expression is:

$$H_{ji}=\max (0,r_{\min }-\Vert x_j-x_i\Vert )$$

(18)

where X_j and X_i are the positions of node j and node i, respectively. The expression of density filtering is as follows:

$${\tilde{\rho }}_j=\frac{1}{{\displaystyle \sum _{i\in N_i}{H}_{ji}}}{\displaystyle \sum _{i\in N_j}{H}_{ji}}{\rho }_i$$

(19)

The unit density ${\rho }_i$ before filtering is called the design variable, and the unit density ${\tilde{\rho }}_j$ after filtering is called the physical density.

It should be noted that when density filtering is used, the objective function and constraint function values in the optimization model are calculated using the physical density, i.e., the filtered density is calculated. When density filtering is applied, the derivatives of the objective function and constraint function concerning the design variables can be obtained by the chain rule.

$$\frac{\partial f}{\partial {\rho }_j}={\displaystyle \sum _{e\in N_j}\frac{\partial f}{\partial {\tilde{\rho }}_i}}\frac{\partial {\tilde{\rho }}_i}{\partial {\rho }_j}={\displaystyle \sum _{e\in N_j}\frac{1}{{\displaystyle \sum _{i\in N_e}{H}_{ei}}}}{H}_{je}\frac{\partial f}{\partial {\tilde{\rho }}_i}$$

(20)

One of the most important aspects in solving the accompanying Equations (14) and (15) is the calculation of the included Zernike coefficients based on the collected surface shape or wavefront data, called Zernike fitting. The circular Zernike polynomial is orthogonal on the unit circle and the coefficients can be solved directly from Equation (21) to obtain the results as follows:

$$\frac{1}{\pi }{\displaystyle \underset{0}{\overset{1}{\int }}{\displaystyle \underset{0}{\overset{2\pi }{\int }}{Z}_{fundi}}}(r,\theta )Z_{fundj}(r,\theta )r\mathrm{d}\theta \mathrm{d}r={\delta }_{ij}$$

(21)

${\delta }_{ij}$ is Kronecker compliance. The coefficient can be solved directly from Equation (22) as follows:

$$Z_{coeff}=\frac{1}{\pi }{\displaystyle \underset{0}{\overset{1}{\int }}{\displaystyle \underset{0}{\overset{2\pi }{\int }}w(r,\theta )Z_{fundi}}}(r,\theta )r\mathrm{d}\theta \mathrm{d}r$$

(22)

where $w(r,\theta )={\displaystyle \sum _{i=1}^N{Z}_{coeffi}{Z}_{fundi}(r,\theta )}$.

It is ensured that the coefficients do not interfere with each other during the calculation process to ensure the accuracy of the calculation results. However, the actual situation is that the surface shape or wavefront information is not a continuous function but discrete data. The surface shape value corresponds to the sampling point $(r_k,{\theta }_k)$ is $w(r_k,{\theta }_k)$. The surface data are expressed as

$$w(r,\theta )={\displaystyle \sum _{i=1}^N{Z}_{coeffi}{Z}_{fundi}(r,\theta )}+\mathrm{\Delta }w(r_k,{\theta }_k)$$

(23)

Supposing the number of sampling points is M. Equation (22) is transformed by $w_k\equiv w(r_k,{\theta }_k)$ as follows:

$$Z_{coeffj}=\frac{1}{M}{\displaystyle \sum _{k=1}^M{w}_k}{Z}_{fundjk}$$

(24)

If the Zernike polynomial maintains orthogonality when the data are discretized, the orthogonality condition in the discrete case becomes as follows:

$$\frac{1}{M}{\displaystyle \sum _{k=1}^M{Z}_{fundik}}{Z}_{fundjk}={\delta }_{ij}$$

(25)

However, the case of limited data sampling points can occur as follows:

$$\frac{1}{M}{\displaystyle \sum _{k=1}^M{Z}_{fundik}{Z}_{fundjk}}\ne \frac{1}{\pi }{\displaystyle \underset{0}{\overset{1}{\int }}{\displaystyle \underset{0}{\overset{2\pi }{\int }}{Z}_{fundi}}}(r,\theta )Z_{fundj}(r,\theta )r\mathrm{d}\theta \mathrm{d}r={\delta }_{ij}$$

(26)

Zernike polynomial after discretization loses orthogonality. Therefore, there is always a theoretical error when calculating the coefficients of the Zernike polynomial using Equation (26). To solve the above problem, the method of Nth order Legendre polynomials [28] is used in this paper to construct the discrete orthogonal Zernike polynomials.

The equation for a polar coordinate $r({\eta }_k)$ is as follows:

$$r_k=\sqrt{\frac{1+{\eta }_k}{2}}\in (0,1),k=1,2,\dots ,N$$

(27)

where ${\eta }_k$ satisfies Equation (28) as follows:

$$P_N({\eta }_k)=0,{\eta }_k\in (-1,1),k\in 1,2,\dots ,N$$

(28)

$P_N$ is a Legendre polynomial of order N. The order N is chosen to be the same as the order of the Zernike polynomial, and Equation (28) shows that ${\eta }_k$ is the root of the equation formed by the Legendre polynomial of order N.

The equation for the polar coordinates $\theta$ is as follows:

$${\theta }_j=\frac{2\pi j}{2N-1},j=0,1,\dots ,2N-2$$

(29)

The coordinates of the discrete points are constructed as follows:

$$X=\{x_{jk}=(r_k,{\theta }_j),k=1,2,\dots ,N,j=0,1,\dots ,2N-2\}$$

(30)

$$\begin{array}{l}\frac{n+\left|m\right|}{2}+\frac{n^{\prime }-\left|m^{\prime }\right|}{2}\le 2N-1\\ \frac{n-\left|m\right|}{2}+\frac{n^{\prime }+\left|m^{\prime }\right|}{2}\le 2N-1\\ n,n^{\prime }\in N\hspace{1em}m,m^{\prime }\in Z\end{array}$$

(31)

When Equation (31) is satisfied, the discrete points have the following properties.

$$\begin{array}{l}{\displaystyle \sum _{k=1}^N{\displaystyle \sum _{j=0}^{2N-2}{Z}_n^m}}(r_k,{\theta }_j)Z_{n^{\prime }}^{m^{\prime }}(r_k,{\theta }_j)\frac{A_k}{2(2N-1)}\\ =\frac{1}{\pi }{\displaystyle \underset{0}{\overset{2\pi }{\int }}{\displaystyle \underset{0}{\overset{1}{\int }}{Z}_n^m}}(r,\theta )Z_{n^{\prime }}^{m^{\prime }}(r_k,{\theta }_j)\rho d\rho d\theta ={\delta }_{n{n}^{\prime }}{\delta }_{m{m}^{\prime }}\\ A_k={\displaystyle {\int }_{-1}^1{l}_k}(x)dx\\ l_k(x)={\displaystyle \prod _{\begin{array}{l}j=1\\ j\ne k\end{array}}^N\frac{x-{\lambda }_j}{{\lambda }_k-{\lambda }_i}}\end{array}$$

(32)

where A_k is the product factor and $l_k(x)$ is the Lagrangian interpolation basis function, ${\lambda }_k={\eta }_k$.

Additionally, by numerical arithmetic example (the surface of a circular optical element with a radius of R, using a 36-term Fringe Zernike polynomial, N = 6), it can be concluded that the errors of the Zernike polynomial coefficients Z_coeffi are of the order of 10–13. By constructing discrete orthogonal Zernike polynomials, the problem of small changes in sampling points leading to large changes in the free-form surface is solved because the Zernike polynomials are not orthogonal in the discrete coordinates inside the unit circle; at the same time, the sampling points required for the surface model are separated from the finite element nodes, which overcomes the problem that the higher the accuracy required for calculating the wave aberration of high-precision optical systems by the finite element method, the more sampling points required and the larger the computational effort.

Equation (2) is approximated as convex programming by MMA, and its specific principle is shown as follows.

$$f_j^n(x)=r_j^n+{\displaystyle \sum _{i=1}^I(\frac{p_{ji}^n}{A_i^n-x_i}+\frac{q_{ij}^n}{x_i-B_i^n})}$$

(33)

$$p_{ji}^n=\{\begin{array}{l}(A_i^n-x_i{)}^2\frac{\partial f_i(x^n)}{\partial x_i},\hspace{1em}\frac{\partial f_i(x^n)}{\partial x_i}>0\\ 0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em},\hspace{1em}\frac{\partial f_i(x^n)}{\partial x_i}\le 0\end{array}$$

(34)

$$q_{ji}^n=\{\begin{array}{l}-(x_i^n-B_i^n{)}^2\frac{\partial f_i(x^n)}{\partial x_i},\hspace{1em}\frac{\partial f_i(x^n)}{\partial x_i}<0\\ 0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em},\hspace{1em}\frac{\partial f_i(x^n)}{\partial x_i}\ge 0\end{array}$$

(35)

$$r_j^n=f_j^{}(x^n)-{\displaystyle \sum _{i=1}^I(\frac{p_{ji}^n}{U_i^n-x_i^n}+\frac{q_{ij}^n}{x_i^n-L_i^n})}$$

(36)

where, B and A are left and right moving asymptotes respectively, p and q are parameters determined by the corresponding B and A, and r is the undetermined parameter. $f_j^n$ is the approximate function of the objective function. By iterating the subproblems of the approximate function, the optimal solution in the tolerance range can be obtained, and this iterative convergence is very fast under reasonable settings.

The flow chart of the optimization model execution procedure is shown in Figure 1.

3. Numerical Examples and Simulation Analysis

The process factor k₁ is the comprehensive influence of illumination, mask, objective, photoresist, exposure method, alignment, and other processes of the lithography system on the resolution [3]. When only the influence of the objective lens on the process factor is considered, it is reflected in the aberration control of the optical lens. With the stringent aberration control requirements, the lithographic projection objective puts limits on static parameters such as lens material uniformity and processing and mounting accuracy. It is also necessary to consider the temporal change characteristics such as air pressure, vibration, and temperature state during the operation of the lithography machine, among which the impact of the thermal effect on the optical path is the most difficult to control.

In the multiexposure technology of the lithography machine, multipole and user-defined lighting modes are adopted. Due to the temperature change in the working environment and the thermal absorption of components caused by multimode light source irradiation, the lens is deformed and the refractive index changes, resulting in the thermal aberration of the projection objective lens changing at any time and many types of modes. Eliminating thermal aberration has become a difficult problem in the development of advanced lithography machines that must be solved.

The size and form of wave aberration caused by different illumination modes are different. A commonly used research method is to decompose the wave aberration of the system into different aberration modes. Under two-pole illumination, the wavefront aberration of the central field of view of the lithography objective is decomposed by Zernike polynomial to obtain the mode of wavefront aberration. In addition to the defocus term, clover aberration and astigmatism are the main components [29].

For simplicity, this paper only takes the astigmatism correction active deformation mirror as an example to study the optimization theoretical method, but the whole method can be extended to the support structure design of any mode wave aberration active deformation mirror. The software platforms used in this manuscript are Matlab, Comsol, and ZYGO interferometer. The initial model and size of the active deformation mirror studied in this paper are shown in Figure 2.

The whole model is divided into two parts: the deformation mirror and its supporting structure. The deformation mirror is a bent moon concave lens, and its structure and main dimensions are shown in Figure 2. To realize astigmatism correction, the supporting structure is composed of an inner ring elastic support structure, an outer ring fixed support structure, and eight flexible plates. The deformation mirror is supported by the inner ring elastic support structure, eight flexible plates connecting the inner ring elastic support structure and the outer ring fixed structure, evenly distributed along the circumference, respectively located in $\theta =0,\pi /4,\pi /2,\hspace{0.17em}\hspace{0.17em}3\pi /4,\hspace{0.17em}\hspace{0.17em}\pi ,\hspace{0.17em}\hspace{0.17em}5\pi /4,\hspace{0.17em}\hspace{0.17em}3\pi /2,\hspace{0.17em}\hspace{0.17em}7\pi /4$, without loss of generality. To correct astigmatism Z4 as an example, the actuators are placed in $\theta =0,\pi$, applying the vertical paper to face the outside force, actuators are placed in $\theta =\pi /2,3\pi /2$, applying the force in the vertical paper, as shown in Figure 2.

Using the discrete points listed in Equation (30) to build a geometric model, the Zernike polynomial surface shape can be calculated more accurately by taking 66 points on the surface of the circular optical element for a 36-term Fringe Zernike polynomial with N taken as six.

The load and boundary conditions of the model are shown in Figure 2. In the figure, the red area is the optimization design area. The inward force perpendicular to the paper surface is applied at boundaries one and two with the size of 117N. The outward force perpendicular to the paper surface is applied at boundaries three and four, with the size of 117N. Boundary five is fixed, and the mirror is used through gravity. The material used in the model is: the mirror body is fused quartz and the flexible support structure in contact with the mirror body is brass. The material properties of the structure are shown in

Table 1. Performance parameters of mirror body and supporting structure materials.

Material	Density (E-6 kg/mm3)	Elastic Modulus (kg f/mm2)	Poisson’s Ratio
Fused silica	2.2	7459	0.167
Brass	8.43	20,000	0.277

. The objective function is the maximum absolute value of Z₄ in the fourth term of the Zernike coefficient. Through the analysis of the initial structure, it can be obtained that Z₄ is negative, and the objective function can be the minimum value of -Z₄, corresponding to Seidel astigmatic aberration.

In order to meet the symmetry requirements of the optimized structure during optimization, the symmetry constraint is applied to make the structure symmetrical along the circumferential direction quarter to quarter. Tetrahedral mesh was adopted with the number of 421,663, and the optimized mesh model is shown in Figure 3.

Since the optimized area is the flexible plate connecting the inner ring elastic support structure and the outer ring fixed structure, to ensure the convergence of the finite element analysis, V_max = 0.7 and V_min are taken as 0.6, 0.5, 0.4, 0.38, 0.34, and 0.3, and $p=0.3$.

Because the optimization model adopts the quarter circumference symmetry model and the optimization results are symmetric, only the quarter circumference symmetry topology optimization results are given here. At the same time, in order to observe the topology optimization results from a better perspective, the fixed structure of the outer ring, the elastic support structure of the inner ring, and the deformable mirror on it are not drawn. The changed results of the topological configuration with volume constraints are shown in the figure. Red indicates the solid structure and blue indicates the materials to be removed.

In this paper, the wave front difference coefficient is taken as the objective function of the topology optimization model, and the constraint conditions are the finite element equation of force balance and volume constraint. The qualitative analysis of the optimization model can be concluded as follows: the smaller the volume constraint value, the smaller the structural stiffness and the smaller the objective function value. However, when the volume constraint is less than 0.1, the continuous structure cannot be obtained. When the volume constraint is greater than 0.5, the continuity structure can be obtained, but the objective function value is large. The volume constraint is between 0.1 and 0.5, and the optimization effect is obvious and general. Therefore, the result of V_min = 0.4 is randomly used to show the optimization process and optimization structure.

When V_min = 0.4, the iterative process is shown in Figure 4.

The iterative process shows that the initial value of the objective function is $-3.85\times {10}^{-4}$, the value of the objective function at the end of optimization is $-1.15\times {10}^{-3}$, and the absolute value of Z₄ is improved by nearly one order of magnitude, so the optimization effect is obvious.

From the topology optimization results, it can be seen that when V_min = 0.6, the inner ring elastic support structure and the outer ring rigid fixed structure are connected by a wedge-shaped flexible hinge at position $\theta =\pi /4$, $\theta =3\pi /4$, $\theta =5\pi /4$, $\theta =7\pi /4$. When V_min = 0.5, the elastic support structure of the inner ring and the rigid fixed structure of the outer ring are at $\theta =\pi /4$, $\theta =3\pi /4$, $\theta =5\pi /4$, $\theta =7\pi /4$, and the wedge angle of the wedge-shaped flexible hinge becomes sharp. When V_min = 0.4, the inner ring elastic support structure and outer ring rigid fixed structure in $\theta =\pi /4$, $\theta =3\pi /4$, $\theta =5\pi /4$, $\theta =7\pi /4$, and the wedge angle becomes zero and the wedge-shaped flexible hinge becomes a rectangular notch flexible hinge. When V_min = 0.38, the rectangular notch flexible hinge becomes a complex flexible hinge. When V_min = 0.34, the depth of the inner flexible hinge in the complex flexible hinge becomes deeper. When V_min = 0.3, the inner flexible hinge in the complex flexible hinge is almost disconnected and a new topology configuration appears. When V_min is less than 0.3, the inner flexible hinge in the complex flexible hinge breaks and the complex flexible hinge becomes a false flexible hinge, which cannot be realized in reality. When V_min is less than 0.2, the finite element analysis results do not converge. The overall analysis shows that When V_min is 0.6 to 03, the flexible hinge changes from a wedge-shaped flexible hinge to a rectangular notch flexible hinge and then to a complex flexible hinge, and the deformation capacity of the structure increases as the minimum volume constraint decreases. The range of corrected astigmatic wave aberration increases with the decrease of minimum volume constraint, but the main topology optimization configuration does not change significantly with the decrease of minimum volume constraint, and no new discontinuities appear. When V_min is 0.3 to 0.2, the inner flexible hinge in the complex flexible hinge breaks and a new topological configuration appears, and the complex flexible hinge becomes a false flexible hinge, which cannot be realized in reality. The objective function decreases with the decrease of the minimum volume constraint, that is, Z₄ increases with the decrease of the minimum volume constraint. The larger the range of corrected astigmatic aberration, the better the optimization effect. Since the thickness h of the deformed mirror is 2R < 0.1 compared to the size of the deformed mirror, the deformation w of the thin plate is also very small compared to the thickness h. To verify the reasonableness of the topology optimization results, the design results were analyzed using the theory of thin plate bending with small deflection. Without loss of generality, a schematic diagram of the forces on the flat deformed mirror microelement in polar coordinates was established, as shown in Figure 5.

The equilibrium diagram of the forces acting on the flat plate in polar coordinates is shown in Figure 5, where the sheer force Q_r, the bending moment M_r, and the torque M_rt at the mirror edge belong to the external forces, the bending moment M_t, and the sheer force Q_t. They need to satisfy the equilibrium equations in three spatial directions, and the actual independent forces are only two, and two of the three are independent elements. According to the theory of bending of thin plate with small deflection, the external and internal forces of the deformation mirror of the plate can be solved to satisfy the following relationship.

$$\begin{array}{l}{M}_r=-D\left[\frac{{\partial }^{2}w}{\partial r^2}+v\left(\frac{1}{r}\frac{\partial w}{\partial r}+\frac{1}{r^2}\frac{{\partial }^{2}w}{\partial {\theta }^2}\right)\right]\\ M_t=-D\left[\frac{1}{r}\frac{\partial w}{\partial r}+\frac{1}{r^2}\frac{{\partial }^{2}w}{\partial {\theta }^2}+v\frac{{\partial }^{2}w}{\partial r^2}\right]\\ M_{rt}=-D(1-v)\left(\frac{1}{r}\frac{{\partial }^{2}w}{\partial r\partial \theta }-\frac{1}{r^2}\frac{\partial w}{\partial \theta }\right)\hspace{1em}\hspace{1em}\\ Q_r=-D\frac{\partial }{\partial r}({\nabla }^{2}w)\\ Q_t=-D\frac{\partial }{r\partial \theta }({\nabla }^{2}w)\end{array}$$

(37)

where $D=\frac{E{h}^3}{12(1-v^2)}\hspace{1em},\hspace{0.17em}w(r,\theta )=Z_5{(\frac{r}{R})}^2\cos (2\theta )$, the equation can be obtained as follows:

$$\begin{array}{l}{M}_r=-\frac{E{h}^3}{6(1-v)}{Z}_5\frac{1}{R^2}\cos (2\theta )\\ M_t=\frac{E{h}^3}{6(1+v)}{Z}_5\frac{1}{R^2}\cos (2\theta )\\ M_{rt}=\frac{E{h}^3}{6(1+v)}{Z}_5\frac{1}{R^2}\sin (2\theta )\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\\ Q_r=0\\ Q_t=0\end{array}$$

(38)

Equation (38) demonstrates that in position $\theta =0,\pi$, the bending moment M_r obtains a great value, the direction makes the deformation vertical paper face outward. At position $\theta =\pi /2,3\pi /2$, the bending moment M_r obtains a minimum value in a direction such that the deformation is perpendicular to the paper facing inward. In position $\theta =\pi /4$, $\theta =3\pi /4$, $\theta =5\pi /4$, $\theta =7\pi /4$, the deformation W($r$,$\theta$) = 0, and M_r = 0. This position corresponds to the position where the deformation is zero, which is the support position. The topology optimization results are consistent with these positions and with the topology configuration designed by Hugot et al. [30] using elastic beam theory.

Considering comprehensively that the astigmatism correction ability of deformation mirror increases with the decrease of minimum volume constraint and the characteristic frequency of flexible hinge decreases with the decrease of minimum volume constraint, we used the topology optimization result of V_min = 0.4 as the reference for the final topology optimization configuration, the main topology configuration was extracted and the structure was designed from scratch, and the results are shown in Figure 6.

The finite element model of the structure designed in Figure 6 was established and divided into tetrahedral mesh with 136,709 meshes and compared with the structure before optimization shown in Figure 2, with the comparison results shown in

Table 2. The comparison results of the structural analysis before and after optimization.

Structure	Z4 (λ = 632.8 nm)	Axial Mode Natural Frequency (Hz)	Z4 Vibration Shape Natural Frequency (Hz)
Structure before optimization	$-1.8\times {10}^{-4}$	891.906529	997.668618
Optimized structure	$-1.3\times {10}^{-3}$	590.523771	656.090888

Then, we used the n-step rotating average method to detect the research object with high precision. Zygo GPI XP/DTM vertical interferometer was used in the surface shape detection experiment, where the working wavelength of the interferometer was 632.8 nm and the resolution of CCD is 1K × 1K. During the experiment, the components were assembled precisely, the temperature control condition was better than 0.01 °C, the vibration condition was better than the VC-E standard, the repeatability of RMS measurement was better than 0.1 nm. During the inspection, the mirror surface shape supported by the structure before and after optimization was detected. The results shown in Figure 7 are the measurements after removing the surface shape errors such as defocus, tilt, coma, and spherical aberration.

From Table 2 and

Table 3. The test results of optimization target Z₄.

	The Test Results of Optimization Target Z4 after Optimization	The Test Results of Optimization Target Z4 before Optimization
Z4 (λ = 632.8 nm)	0.001	0.008

and Figure 7, it can be concluded that the absolute value of Z₄ of the structure designed based on the topology optimization results was improved by nearly one order of magnitude compared with the absolute value of Z₄ of the structure before optimization, At the same time, the vibration characteristics of the optimized structure meet the design requirements, which verifies the effectiveness of the optimization model, and the vibration characteristics of the optimized structure meet the design requirements, the optimization effect is remarkable.

4. Conclusions

Firstly, according to the demand for a high-precision wavefront compensation to realize the whole field of view aberration compensation process using an active free-form surface, the objective function and design constraints of the optimization model were constructed based on the Zernike coefficient. On this basis, the wavefront aberration topology optimization model with tens of thousands of design variables was established using the topology optimization method of the variable density method. The sensitivity of Zernike coefficients to topology-optimized design variables is derived using the concomitant variables method, which overcomes the computationally intensive problem of solving the sensitivity in the topology optimization process. To address the problem that Zernike polynomials do not have orthogonality in discrete coordinates inside the unit circle, and small changes in sampling points lead to large changes in the solution of the free-form surface and the accompanying equations, discrete orthogonal Zernike polynomials were used to characterize the free-form surface of the actively deforming mirror, which ensured the computational accuracy and simplified the computational process by realizing the characterization of the free-form surface of the actively deforming mirror with small samples.

Then taking a light wave aberration correction of the deformable mirror machine system as an example, under the theoretical framework of the finite element numerical discrete, we adopted a finite element basis function unit and a numerical integral program implementation solution of structure deformation and Zernike coefficients in weight constraint conditions, based on a mobile asymptote algorithm, the optimization model into the solution. The optimal topology configuration of the deformation mirror satisfying the convergence condition is obtained. Using thin plate bending theory with small deflection, the rationality of topology optimization results was analyzed from the perspective of theoretical analysis. At the same time, the topology optimization results are consistent with the topology configuration designed by Hugot et al. using elastic beam theory, which verifies the effectiveness of the algorithm.

Finally, the topology optimization results at V_min = 0.4 were used as a topology reference to optimize the final configuration, extract the main topology configuration, and analyzed the optimization effect and characteristic frequency using the finite element analysis, and the mirror surface shape supported by the structure before and after optimization was detected. The following results were obtained from the simulation results and experimental results. The absolute value of the objective function Z₄ was improved by nearly one order of magnitude compared with the absolute value of the structure Z₄ before optimization, while the vibration characteristics of the optimized structure met the design requirements and the optimization effect was remarkable. The optical performance of the deformed mirror was improved, which provided a new scheme for the design of the deformed mirror structure and had a certain engineering practical value.

The algorithm in this paper can optimize the general free-form surface structure topology optimization model with Zernike coefficients constrained by linear combination and has a certain generality. Meanwhile, the Zernike coefficient of the wavefront in image space can be considered as a target to optimize the support structure of the optical system, so as to realize the topology optimization at the system level.