A hybrid eigenmode restoration algorithm for computational lithography problems based on mode matching principle

The mode matching method is an efficient solution to conventional microwave waveguide problems for its dimensionality advantages. It is theoretically extendable to computational lithography problems in nano scales with periodic boundary condition. The high computational complexity of conventional mode matching limits its application significance in highly complicated nano scale problems. This paper introduces an optimized efficient mode matching method for computational lithography problems whose complexity is O(N1.5). This bottleneck is breached by transforming governing equations with Lanczos algorithm. A novel hybrid eigenmode restoration algorithm is proposed to solve the eigenmode accuracy loss by involving the Arnoldi method. Benchmark cases verify the accuracy of the restored eigenmode and mode matching method with from microwave to nano scale problems. The efficiency is further optimized by exploring the relevance between iteration step and structure complexity. The non-uniform iteration step is introduced for efficiency optimization. A real mask component is modelled as a multi-junctional waveguide structure with periodic boundary condition. Simulation results validate the effectiveness of the proposed method and efficiency benefits are discussed through comparison with other solvers including HFSS and RCWA. Results indicate that the proposed method possesses good flexibility and application significance for computational lithography problems with high complexity.


I. INTRODUCTION
In recent years, there are growing requirements and achievements in numerical techniques for inhomogeneous dielectric waveguide analysis. These techniques are also extendable from microwave to nanoscale computational lithography research for further enabling Moore's law in next a few years [1][2][3][4][5] . The mode matching method is one of the effective solutions for these problems [6,7] . It evolves from Ez-Hz formulation to transverse electromagnetic field formulations to solve spurious mode problems. Unknowns are only defined on transverse interfaces at the structure discontinuity. This feature overcomes the curse of dimensionality by transforming a large three dimensional problem into several two dimensional or 2.5 dimensional problems [8,9] . Numerical modelling techniques enable mode matching methods to solve complicated waveguide structures using finite element or finite difference techniques [10,11] . Existing solvers are capable of solving three dimensional on-chip interconnects with two dimensional time and memory consumption [12][13][14] . The dimensionality advantage makes mode matching method extendable to nanoscale problems. Transverse field definitions are also suitable for modelling multi-layered structures, which are common in computational lithography analysis. Existing works verify that the periodic boundary condition (PBC) enable mode matching method simulate periodic structures [15] . For nonperiodic structures, the mode matching method is also applicable with a sufficiently large buffer zone and PBC. These characters provide flexibility in analyzing various computational lithography problems.
In mode matching, electromagnetic fields are expanded by a series of eigenmodes, which are calculated and matched at each junction discontinuity. The kernel of the mode matching method is the eigenmode computation and expansion. Existing mode matching methods usually acquire eigenmodes by spectral decomposition techniques. Its computational complexity is usually O(N 3 ) [16] in which N is the unknown number for transverse electromagnetic fields. In many computational lithography problems, the high geometry complexity and nanoscale wavelength results in a large N. In these cases the O(N 3 ) complexity becomes an efficiency bottleneck in mode matching implementation. This gets more prominent for multilayered structures. An implicit mode matching (IMM) algorithm based on Lanczos method [17,18] is proposed to solve this efficiency limitation. It reformulates governing transverse wave equations into implicit form to avoid the explicit eigenmode solving procedure. According to the boundary condition, the matching procedure is applied on transverse field directly instead of eigenmode expansion coefficients. The IMM method could reduce the complexity from O(N 3 ) to O(N 1.5 ). However, the price of efficiency enhancement is the loss of accuracy. The non-strict orthogonalization procedure in Lanczos iteration leads to the eigenmode accuracy deviation, which gets non-negligible for highly complicated structures. Another limitation is its high dimensional block matrix in solving multilayered structures. Moreover, the IMM makes it difficult to interpret eigenmode behaviors such as propagation, reflection and decaying effects. This also constrains its physical essence exploration. Therefore, even though efficiency benefits are distinctive, the application significance of IMM in computational lithography problems is still unclear.
Works presented in this paper are motivated by developing an efficient mode matching algorithm aiming at analyzing complicated lithography problems. Constrains from IMM accuracy are solved but preserving its efficiency advantage. The eigenmode accuracy problem is corrected by proposing a new hybrid eigenmode restoration algorithm based on Arnoldi method [19,20] . Governing equations are also transformed into explicit form using restored eigenmodes. Accuracy of eigenmode restoration and new mode matching procedures are validated using benchmark cases with both perfect electric conductor (PEC) boundary condition and PBC. The High-Frequency Structural Simulator (HFSS) [21,22] software and RCWA (Rigorous Coupled Wave Analysis) method [23][24][25][26] are adopted to verify the proposed method. Relevance between Lanczos iteration steps and structure complexity are explored by specific simulation cases. The analysis results in a non-uniform iteration strategy to further elevate efficiency. A mask component with high geometry complexity is designed according to practical industrial analysis examples. Simulation results and comparisons with HFSS, RCWA and IMM validate the performance in both accuracy and efficiency. These consolidate the proposed method in solving complicated computational lithography problems. This paper is organized as follows. Section II presents governing equations and the hybrid eigenmode restoration strategy. The mask structure modelling using the proposed mode matching method is introduced in section III. Results and discussions in section IV verify the effectiveness of the proposed method from various aspects.

A. PROBLEM FORUMULATION
Governing wave equations for transverse electromagnetic fields under the finite difference discretization are [27] : In inhomogeneous medium, Le and Lh are sparse and asymmetric operators. Two operators are approximated in matrix form with finite difference modelling. As the Z component of electric and magnetic fields could be derived from transverse components according to Maxwell equation, transverse electric and magnetic fields are taken as unknowns. In following contents vectors E and H are taken to represent transverse electromagnetic fields, mediums are assumed to be isotropic and non-dispersive so that tensorsˆs μ and s ε  are simplified to scalars  and  .
In conventional mode matching, eigenvalue problems in (3) are solved using numerical techniques like spectral decomposition with O(N 3 ) complexity. This is the origin of its efficiency limitation. In contrast, equations (3) are solved in an alternative way in IMM by introducing the Lanczos method. It is based on Krylov subspace method [28,29] to find basis vectors to approximate the original high dimensional matrix. For the original matrix Le, the Lanczos method generates basis vectors by making an initial guess on the first basis vector and following the iteration procedure: The procedure in (4) starts from an initial guess of v0, the new basis vector is generated by: In Lanczos method it is assumed that basis vectors possess the orthogonality property. Therefore, it is deduced that: Procedures in (4) are transformed into another form: Above iteration procedures are presented in a matrix form: The term indicates a numerical error. Each column of matrix Vk represents a basis vector v. Equation (8) is equivalent to a numerical approximation of Le matrix by a lower dimensional basis space Vk expanded by coefficients in the tridiagonal matrix. Above procedures are major steps for Lanczos algorithm. The Bi-Lanczos method is the extension of Lanczos to solve asymmetric matrix problem by introducing an extra iteration, its implementation is similar to equations (4) to (7). The tridiagonal matrix could be denoted as matrix T with dimension M, which is usually scaled as c· N 0.5 according to Lanczos method's principle, in which c is a constant [30] . The dimension relationship among different matrices could be intuitively interpreted in Figure 1. The matrix D represents the numerical error, which is controllable with increasing iteration steps. The original Le operator could be approximated in an alternative form: Columns vectors in matrices V and W are Bi-Lanczos basis vectors. The dimension of matrices V and W is N-by-M. The numerical term D is neglected. Based on the assumption that most matrix functions adopted in engineering problems could be approximated by a high order Taylor series, equation (9) evolves into a generalized form like: In which f(·) indicates a function of matrix. The multiplication of a matrix function and a vector is the kernel computation component in the mode matching method. It is also the most time-consuming part in the overall procedure. Conventionally the computation of a matrix function costs massive resources. The introduction of the Lanczos algorithm simplifies this part with an alternative approach. For a matrix-vector multiplication in the form of   e f  L v , its computation starts from selecting vector v as an initial basis vector. By following Lanczos iterative procedures, matrix V and T are generated and the multiplication is transformed into: The vector e1 extracts the first column of a matrix. The matrix T is decomposed by spectral decomposition method: The matrix Λ represents a diagonal matrix containing The basis vector qi from the i th column of matrix Q represents an eigenvector corresponding to i  . Equation (11) could be transformed into following form by introducing (12): In ( Procedures from (9) to (13) consumes most resources, they determine the complexity of the mode matching method. The spectral decomposition of tridiagonal matrix T consumes most computational resources with complexity of O(M 3 ). According to the relationship of M~c·N 0.5 , the complexity is equivalent to O(N 1.5 ). This is the computational complexity of IMM. Therefore, the Lanczos algorithm reduces the computational complexity from O(N 3 ) to O(N 1.5 ). This elevates the simulation efficiency for problems with large unknowns. Procedures from equation (4) to (8) are also applicable to Lh and H. As Lanczos and Bi-Lanczos methods are based on the same principle with slight difference, to avoid confusion, only the term Lanczos method is adopted in following contents.

B. HYBRID EIGENMODE RESTORATION
According to the relationship between Le and T, diagonal elements of matrix Λ well approximate eigenvalues of Le.
propagating modes if 2 0 z k  . However, the accuracy of this approximation is limited to the first a few dominant modes due to the orthogonality loss in Lanczos iteration. This accuracy loss is not obvious for simple structures as first a few dominant eigenmodes could well interpret mode physics and mode matching accuracy is not prominently influenced. However, in the case of high geometry complexity, higher order modes are required to capture complicated mode physics and the eigenmode accuracy loss usually results in diverged numerical errors. Therefore, physically i  V q is a reasonable eigenmode approximation but numerically it could not meet accuracy requirement for complicated geometries. This limits its capability in complicated lithography problem analysis.
Above analysis indicates that the source of eigenmode approximation error is the non-strict orthogonality iteration strategy in Lanczos. However, this iteration strategy is also the reason that the problem of   e f  L v could be solved with O(N 1.5 ) complexity. This seems a contradictory in Lanczos algorithm between accuracy and efficiency. Works presented in this paper are motivated by seeking an approach which provides more accurate eigenmode approximation while maintaining Lanczos efficiency benefit. The solution is introducing the Arnoldi method into eigenmode construction. Arnoldi and Lanczos methods are all derived from the Krylov subspace theory [31,32] to project the original high dimensional matrix into a lower dimensional space. However, their basis vector generation procedures are different. In Arnoldi, the iterative procedure for basis generation is: Compared with (4), the distinctive difference of two methods is about the second term of the right-hand side. In Arnoldi a modified Gram-Schmidt orthogonalization procedure is adopted. The e j  L u term is orthogonalized against all previously generated basis vectors from u0 to uj. In contrast, Lanczos method only considers neighboring basis vectors in orthogonalization. This is the reason that Lanczos iteration is a non-strict orthogonalization procedure, which results in basis orthogonality loss. The Arnoldi method does not have this orthogonality loss, but it consumes more resources as all previously generated basis vectors are considered in new basis vector generation. The new basis generation in (5) and the property in (6) are also applicable in Arnoldi method. The difference is about the approximation matrix. Based on (14), a new coefficient is defined as: According to basis orthogonality property of u: Therefore, the matrix Le could be approximated by basis vector u in a similar form like (8). According to (15), the coefficient matrix is not tridiagonal but in a upper Hessenberg matrix form: 11 12 1 Each column of Uk is an Arnoldi basis vector, and all basis vectors possess strict orthogonality. The symbol Θ indicates a submatrix with all zero elements.
Our contribution in hybrid eigenmode restoration is derived from the further exploration of physical essences for matrices V and U. As matrix T from Lanczos method contains dominant eigenspace information, its eigenvalues and eigenvectors are preserved. However, the projection matrix V is replaced by U. The reason of this replacement is based on their similar physical essence, which is projecting the original matrix into a specific Krylov subspace. Therefore, this replacement is for approximating the original problem space using an alternative subspace projection strategy, and i  U q is the restored eigenmodes adopted in following mode matching procedures. As matrices U are generated and stored in an off-line manner, the strict orthogonalization procedure has limited adverse impact on its efficiency. The adoption of U makes basis i  U q preserve consistent orthogonality in subspace spanning, and higher order eigenmodes could be interpreted more accurately.
As larger eigenvalues are related to more dominant eigenmodes, matrix Λ is sorted in descending order in numerical implementation. The numerical error of the restored eigenmode is quantified as: The operator 2  indicates a L-2 norm computation of the vector. If the error   i  is larger than a specific threshold, the corresponding eigenmode is not selected in mode matching. This usually occurs for higher order eigenmodes with relatively smaller eigenvalues. In this paper the error threshold value is set to 10 -4 . The number of selected eigenmodes M' is related to the energy of the system, which is defined as the squared sum of eigenvalues in matrix Λ. In practical implementation, an energy threshold is specified to define the parameter M'. For example, an energy threshold of 90% indicates that first M' eigenmodes whose eigenvalues in Λ takes up 90% energy of the system. Eigenmodes with   i  larger than the error threshold are removed. For simple structures, first a few dominant eigenmodes could interpret the physical essence of the system, the energy threshold could be relatively small like 70% or 80%. For highly complicated structures, higher order eigenmode contributions are non-negligible and more eigenmodes should be selected in mode matching. In this case, energy threshold value should be large and sometimes the value is VOLUME XX, 2017 9 set to 100%. More discussions about this topic are given in section IV with specific problems.

C. BOUNDARY CONDITION
An eigenvalue problem could not have unique solutions until boundary conditions are applied. Most conventional mode matching methods adopt the PEC boundary condition to confine problems for waveguide structures. PEC is also applicable to the proposed method to solve classic waveguide problems. However, for regular computational lithography problems the PEC boundary condition is not adopted due to improper physical meaning interpretation. The PBC [33][34][35][36] is taken for lithography problems. A two-dimensional PBC is presented in Figure 2 with mathematical presentation: Parameters a and b indicate periods in X and Y directions. Symbols kx and ky are wave numbers. The non-italic symbol j is used to denote the imaginary number whose square is -1. For many lithography problems like mask component analysis, their corresponding geometry configuration might not be periodic. However, the PBC is still applicable by defining a sufficiently large period as the buffer zone. Another prominent benefit of PBC is that the corresponding Floquet modes naturally approximate the plane wave excitation. This provides a more reasonable description for most computational lithography problems.

III. MASK STRUCTURE MODELLING USING MODE MATCHING METHOD
In this paper the mask component is selected to verify the capability of the proposed method. A mask without substrate could be modelled as a two-junctional waveguide structure with PBC, as illustrated in Figure 3. Junctions I and III are free space. The blue area in section II represents a dielectric region of component. The height of the dielectric region in Z direction is usually much smaller than its size in XOY plane. For structures with a substrate, or more complicated multilayer structures, the proposed method is still applicable by modelling the problem as a multi-junctional waveguide structure. Unknown definitions are presented in Figure 4 for a generalized waveguide structure. According to the Maxwell's equations, transverse components of magnetic field could be derived from transverse electric field from (1) and (2): By utilizing the duality principle, the transverse electric field are presented in a similar manner: Equations (20) and (21) could be presented in a compact form by introducing operators A and B: Le and Lh operators could be depicted with more compact form using A and B by utilizing some mathematical tricks:  Both transverse electric and magnetic fields are matched in a similar manner at the discontinuity due to the duality principle. In following contents only electric field matching equations are presented. According to Figure 4, the electric field matching at the first junction results in: The matrix Mm,n is defined as: Equation (26) is expanded using restored eigenmodes as: The matching in middle junctions takes a similar procedure: Equation (29) is transformed to the explicit mode matching forms by expanding electric fields with restored eigenmodes: The transmitted fields at the last junction are matched as: The elimination of transmitted field results in: The explicit eigenmode expansion transforms (32) into: A block matrix is constructed by reorganizing (28), (30) and (33). Equation (34) The first row indicates the matching at the first junction interface. The element 11 Ψ indicates a matrix corresponding to the first left hand side term of (28). The symbol Θ indicates a sub-matrix with all zero elements. The symbol 0 on the right-hand side indicates a vector with all zero elements. Rows 2 to 5 correspond to the middle junction matching. The last row denotes the mode matching at the last junction. Unknown dimensions of the matrix system in equation (34) are prominently reduced compared with IMM. Moreover, the block matrix is well conditioned which assists the iterative solver converge fast to calculate expansion coefficients. The overall computational complexity is O(N 1.5 ). Since eigenmodes need to be stored in the memory for block matrix construction, the overall storage complexity is higher than but the consumption is still in reasonable scale.

IV. RESULTS AND DISCUSSION
In this section, numerical results of many simulation cases are discussed to verify the proposed method from various aspects. The first subsection validates the accuracy of the hybrid eigenmodes using representative waveguide structures. Both microwave and nano-scale structures are adopted in the second subsection to validate the mode matching method accuracy and its capability of processing arbitrarily shaped structures. The relevance between Lanczos iteration steps and structure complexity is explored and discussed in the third subsection to support the non-uniform iteration strategy. A mask structure with high complexity is presented in the fourth subsection. Results from various solvers are presented to evaluate the accuracy and efficiency of the proposed method. Numerical complexities of the proposed method and RCWA are discussed in the last subsection to prove the benefits of the proposed method for specific problems.

A. EIGENMODE ACCURACY VERIFICATION
A dielectric rectangular waveguide is selected to verify restored eigenmodes. The waveguide possesses an aspect ratio of 2 and relative permittivity of 2.25 in the center area. The PEC boundary condition is applied. Electric field intensity distribution of eigenmodes Ex11, Ex21, Ex12 and Ex22 are presented in Figure 5. Dotted lines mark the boundary between dielectric core and the cladding. Good matches are presented compared with eigenmodes calculated from the analytical solver in [7]. Normalized propagation constant Ps and the normalized frequency B are defined according to [7] to further verify the eigenmode accuracy: The normalized frequency B for eigenmodes in Figure 5 are all 2.5, and the corresponding normalized propagation constant Ps are 0.83,0.74,0.46 and 0.37 respectively. Figure 6 present dispersion curves for aspect ratios of 2 and 1.5 with a medium contrast of 2.25. Good matches could be observed. The root mean square error (RMSE) are 0.007 and 0.01 respectively compared with results in [7].

B. MODE MATCHING METHOD VERIFICATION
As both PEC and PBC boundary conditions are applicable on the proposed method, the mode matching verification in this section is divided into two parts. The first part takes the PEC boundary condition to model a microwave waveguide structure with a complicated geometry pattern. The second part uses PBC to analyze a periodic structure.
The geometry configuration of an arbitrarily shaped waveguide is presented in Figure 7. All length units are inches. This model is designed according to [37]. A cross iris is embedded between two circular waveguide cavities. Each cavity is also connected with a rectangular waveguide as input and output ports. This structure is modelled as a four junctional waveguide component. Results from HFSS software are taken as verification. Figure 8 present magnitude and phase of S11 and S21 parameters with eigenmode energy threshold of 100%. Good matches verify the capability of the proposed method in solving highly arbitrarily shaped waveguide structures. Numerical deviations compared with HFSS are evaluated by using normalized RMSE for magnitude and phase of S11 and S21. Contents in Table I are normalized RMSE with eigenmode energy threshold of 40%, 60%, 80% and 100%. As higher order eigenmodes are usually excited by complicated structures, larger energy threshold is selected to guarantee the accuracy.  The profile of the periodic structure in transverse dimension is presented in Figure 9, all length units are um. The structure is periodically deployed in horizontal and vertical dimensions. Periods are 1um and 0.5um as marked in Figure 9. The thickness of the structure is 0.1um. Relative permittivity and permeability are 10 and 1 respectively. The planewave is incident from normal direction. Wavelengths range from 1um to 2um. Considering the simplicity of the structure geometry pattern, the energy threshold value in eigenmode filtering is set to 70% to meet the 10 -4 error tolerance. This is beneficial for efficiency enhancement while maintaining the accuracy. Reflection and transmission coefficients are presented in Figure 10 with TE and TM polarizations. The TE polarization is defined for the case in which incident electric field is along the horizontal axis (The direction whose length is 1um). The TM polarization is correspondingly defined for the electric field along vertical axis. Results are compared with RCWA solver and HFSS. Good results consistency is presented. The RMSE of coefficients compared with HFSS is defined to quantify the numerical accuracy of IMM and the proposed method. Table  II presents RMSE with the same iteration steps. Both methods present higher accuracy with larger iteration steps. Under the same iteration step, the proposed method could reflect smaller RMSE due to the basis orthogonality correction by introducing Arnoldi projection matrix.

C. LANCZOS ITERATION STEP DISCUSSION
The Lanczos iteration step number M determines the overall efficiency as well as eigenmode accuracy. Theoretically M is scaled as c N  . However, it is difficult to evaluate the specific constant c for a particular problem. Even though c does not influence the scale relationship between M and N, its impact on efficiency is non-negligible. Existing analysis indicate potential relevance between c and structure complexity in geometry and medium contrast. This section further explores this relevance.
Theoretically a larger Lanczos iteration number provides better approximation in (13) with more accurate eigenmodes. However, a quantitative descriptor is required to interpret how large the iteration steps should be. The self-convergence iteration step is introduced as a complement for numerical performance evaluation. It is defined as an iteration step that could make the residual error between left and right hand sides of (9) smaller than a specific tolerance. Connections between geometry complexity and iteration step are studied by designing three structures with increasing geometry complexity as in Figure 11. The black region is filled with homogeneous medium. Unknown numbers are all 5,000 and the residual error tolerance is 10 -4 . Table III presents the selfconvergence iteration step of three structures with different medium contrasts. Contents in the table are the statistical average values. Both geometry complexity and medium contrast present positive correlation relevance with iteration steps. Therefore, when the structure possesses high complexity in geometry or medium configuration, more Lanczos iterations are required for higher order eigenmode generation to capture the physical essence. This iteration step character could be explained in an alternative and more intuitive manner. By generating the system matrix Le matrix for three structures, matrix condition numbers could be analyzed. For the medium contrast of 10, condition numbers for three geometry patterns are 1.74x10 4 , 8.57x10 4 and 3.65x10 5 respectively. It is well known that the larger condition number usually indicates the dominance of larger eigenvalues. As Lanczos method is based on Krylov subspace method, it usually requires more iteration steps for new subspace spanning to overcome the dominance of the large eigenvalue. Therefore, the relevance between Lanczos iteration number and geometry complexity is similar to the mechanism that iterative solvers usually need more iteration steps to solve matrix problems with larger condition numbers. According to above analysis, the uniform Lanczos iteration step adopted by conventional mode matching method is not reasonable, especially for multi-layer structures with different geometry complexity patterns. For junctions with simple geometry configurations, the excessive iteration step is a waste of computational resources, while insufficient iterations could hardly satisfy accuracy requirement for junctions with high complexity. A more intelligent iteration strategy is favorable for efficiency optimization. As selfconvergence iteration step could be evaluated with a specified residual error tolerance, an alternative iteration strategy is preconfigured before simulation. For a problem like Figure 3, the first junction is composed of free space and its iteration step is much smaller than that of the second one. The structure in Figure 11(c) is taken to verify this strategy. Iteration step configuration and computation time for two strategies are presented in Table IV. The non-uniform iteration strategy saves over 50% of time with minor impact on accuracy. For multi-layered structures benefits of the nonuniform iteration strategy get more prominent.

D. SIMULATION AND DISCUSSION FOR A COMPLICATED MASK COMPONENT
Previous discussions and results indicate the potential of the proposed method in solving computational lithography problems. In this subsection a mask with high geometry complexity is introduced for verification. Figure 12 presents a mask component from a practical industrial designing example. The right zone of the figure amplifies the part confined by red dotted lines in the left one. Some geometrical parameters are marked. It could be noticed that layout in Figure 12 is periodic, which makes the adoption of PBC reasonable without the necessity of introducing the buffer zone. According to Figure 12, a mask structure is extracted with slight geometry deviation as presented in Figure 13. The black region is composed of Molybdenum Silicide medium whose relative dielectric constant is (2.43+j · 0.6) 2 . Length, width and thickness parameters are 2245nm, 708nm and 68nm. The unknown number is 215,000 with 193nm wavelength. The structure is excited by normally incident plane waves with TE polarization, where the electric field is along the direction of horizontal axis in Figure 13. The non-uniform iteration strategy is adopted. Junction II takes 3,000 iteration steps. Junctions I and III with free space only needs 200 iteration steps. The transmitted electric field distribution in the lithography region is presented in Figure  14. The field intensity distribution presents the consistency with mask geometry pattern. Transmission coefficients and the phase of total transmitted electric field are presented in Figure 15 with different wavelengths. Results are also compared with IMM and RCWA methods. Good matches could be visually observed to verify the simulation accuracy. The procedure is also extendable to mask component with multi-layer substrate structures by modelling as multiple junctional waveguide problems.  Compared with similar accuracy performance, difference in efficiency among methods are more prominent. The RCWA method is based on Bloch-Floquet theorem and naturally suitable for solving problems with periodic property.
It is mathematically elegant by transforming the problem into the Fourier domain. However, RCWA is intrinsically still a mode matching method by matching transverse fields at the junction discontinuity. Fields are still represented in eigenmode expansion form and the only difference is that eigenmodes are expanded in Fourier domain, which are still solved in RCWA with O(N 3 ) complexity. In contrast, IMM and the proposed method possess a computational complexity of O(N 1.5 ).
Contents in Table V and VI are time and memory consumption for three methods under the wavelength of 193nm and 217nm. The self-convergence iteration step for Lanczos is 3000 and RCWA method needs 3721 modes. Large time consumption of RCWA is caused by its high complexity in eigenmode solving. The proposed method reduces computational time contributed by block matrix dimension reduction and non-uniform iteration strategy. As eigenmodes are explicitly stored in the memory for the proposed method, the corresponding memory consumption is higher than IMM.

E. NUMERICAL COMPLEXITY
There is the necessity to give extra discussions between RCWA and the proposed method in numerical complexity. Even though RCWA is widely used in many periodic structure analysis, its capability in generalization to more complicated structures is still not clear. Numerical studies indicate that RCWA method is not a suitable option for two kinds of problems. The first one is related to geometry complexity. For simple structure like a rectangle, the RCWA could use small amount of plane waves (eigenmodes) to approximate transverse fields with good accuracy. However, much more plane waves are needed to capture physics of detailed structure characters with high complexity. Due to its O(N 3 ) complexity in solving eigenmodes and the quadratic mode number increment character, its resource consumption in time and memory does not reveal any advantages compared with other methods. The second adverse case is structures composed of multiple junctions. For a structure with M1 junctions, the complexity of RCWA is   for solving the block matrix problem in equation (34) and N2 is the unknown number with finite difference modelling. This means RCWA method needs to conduct   3 1 O N computation for M1 times. Even though M2 is apparently larger than M1, the overall complexity is maintained at O(N 1.5 ) level. Therefore, for problems with high complexity in XOY and Z dimensions, the proposed method still presents higher degree of flexibility to balance efficiency and accuracy performance.

IV. CONCLUSION
This paper introduces a novel hybrid eigenmode restoration and explicit mode matching method to solve computational lithography problems. The method is based on IMM but formulations are transformed to explicit mode matching forms. The eigenmode accuracy problem caused by basis vector orthogonality loss is solved by a hybrid restoration algorithm through combination of Lanczos and Arnoldi methods. The overall computational complexity is still maintained at O(N 1.5 ) level. Restored eigenmodes and mode matching method are verified through many simulation cases. The relevance between Lanczos iteration step and structure complexity is explored and a non-uniform Lanczos iteration strategy is proposed for efficiency optimization. Results for a highly complicated mask structure are presented with good accuracy and efficiency performance. Further discussions on computational complexity verify the application potential of the proposed method for highly complicated lithography problems. However, it is necessary to address that the current method still has limitations. For smaller scale lithography problems like 3nm craft which is available right now, governing equations adopted in the proposed method could not interpret the lithography problem comprehensively. At 3nm scale or smaller, more physical phenomena like quantum effects are reflected and the computational lithography problem has evolved into a multi-physics problem. Future studies on modelling smaller scale problems with multiphysics principles might elevate the significance of eigenmode technologies to serve lithography industry.