Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group
Journal Home About Issues in Progress Current Issue All Issues Feature Issues

Three-dimensional plasmonic lithography imaging modeling based on the RCWA algorithm for computational lithography

Huwen Ding, Taian Fan, Libin Zhang, Yayi Wei, and Tianchun Ye

Open Access Open Access

Abstract

This paper reminds the principle and characteristics of plasmonic lithography, and points out the importance of establishing a fast and high precision plasmonic lithography imaging model and developing computational lithography. According to the characteristics of plasmonic lithography, the rigorous coupled-wave analysis (RCWA) algorithm is a very suitable alternative algorithm. In this paper, a three-dimensional plasmonic lithography model based on RCWA algorithm is established for computational lithography requirements. This model improves the existing RCWA algorithm, that is, deduces the formula for calculating the light field inside the structure and proposes the integration, storage and invocation of the scattering matrix to improve the computation speed. Finally, the results are compared with commercial software for the two typical patterns. The results show that the two calculation results are very close, with the root mean square error (RMSE) less than 0.04 (V/m)2. In addition, the calculation speed can be increased by more than 2 times in the first calculation, and by about 8 times by integrating, storing and invoking the scattering matrix, which creates conditions for the development of plasmonic computational lithography.

© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

In traditional lithography techniques, the resolution is usually improved by continuously reducing the exposure wavelength λ and increasing the numerical aperture NA. However, light sources with shorter wavelengths and lenses with larger numerical apertures require significant investment. And because of the diffraction characteristics of light, traditional optical lithography faces the diffraction limit of resolution (∼ λ/4) [1]. As an alternative solution, plasmonic lithography can break the limit of optical diffraction in principle and can be compatible with traditional photolithographic materials and processes without introducing complex and expensive optical lens and short wavelength source. It has gradually developed into a new nano-optical processing technology with high resolution and low cost [2]. Plasmonic lithography is also known as surface plasmon lithography or surface plasmon polariton lithography. The evanescent waves at the mask are resonantly amplified and participate in imaging by exciting surface plasmon polariton (SPP) at the metal/dielectric interface. Finally, exposure and imaging are performed in photoresist to achieve pattern transfer [3–5]. Due to the use of evanescent near-field imaging containing high-frequency information, plasmonic lithography can break the diffraction limit in traditional lithography. Experiment results demonstrated that even with a source wavelength of 365 nm, the imaging resolution can reach about 22 nm (∼1/17 light wavelength) under the condition of single exposure, and it is possible to have further improvement [6,7].

Figure 1(a) shows a typical plasmonic lithography imaging system, which consists of mask substrate Glass, Cr mask, spacing layer PMMA, metal Ag layer, photoresist PR layer, and reflective Ag layer in sequence along the direction of light wave incidence. It can be seen that this imaging system is very different from traditional DUV and EUV lithography systems, as shown in Fig. 1(b). The imaging system is stacked with nanoscale thickness films rather than optical lenses, so the entire imaging area is in the near-field range and the evanescent wave containing the high frequency objective information takes part in imaging to improve imaging quality and lithographic resolution.

 figure: Fig. 1.

Fig. 1. Schematic diagram of imaging principle of plasmonic lithography and conventional projection optical lithography: (a) plasmonic lithography, (b) DUV lithography.

Download Full Size | PDF

To quantitatively and systematically describe the imaging process, obtain the imaging results and further guide how to optimize the process parameters and improve the imaging performance, it is necessary to establish the imaging model of plasmonic lithography. The research of plasmonic lithography imaging model can be divided into numerical method, analytical method and machine learning method. The numerical methods include FDTD [8,9] and FEM [10–12], the analytical method is mainly based on the optical transfer function (OTF) method [13–15], and the other is machine learning method [16]. Each of these methods has advantages and disadvantages, and it is necessary to choose the appropriate method according to different situations and needs. The numerical method is relatively accurate, and many of them have developed commercial software. However, the calculation speed is slow, especially for three-dimensional models. The analytical method represented by OTF method can intuitively explain the physical reasons behind the results and the calculation speed is fast. However, due to the approximation of the mask, the accuracy of imaging results has decreased. Compared with the rigorous numerical method, the fast model based on machine learning and other methods improves the calculation speed at the expense of a certain accuracy. However, the inherent inexplicability of machine learning and other methods generalizes and represents the physical process in the imaging process in a “black box” way, which makes hard to study the physical process in depth, so it is difficult to explore the imaging mechanism through these methods.

As a semi-analytical and semi-numerical method, the rigorous coupled-wave analysis (RCWA) method, also known as Fourier mode method (FMM), is suitable for the simulation of multi-layer structures with transverse non-uniform and periodic patterns incident by plane wave sources [17–19]. And the typical structure is the plasmonic lithography structure shown in Fig. 1(a). In this method, the distribution of the electromagnetic field and the material parameters of the multilayer structure are expressed as the Fourier expansion of infinite series, but the determined truncation order is often needed as the approximate solution, so the speed and accuracy of this method depend on the number of truncation orders. The RCWA method takes into account the advantages of numerical method and analytical method, which not only ensures high calculation accuracy, but also has fast running speed. In addition, the method is based on the solution of rigorous Maxwell equations, conforms to the physical principles and laws, and avoids the inexplicability of “black box”.

Since it was proposed by M. G. Moharam and T. K. Gaylord in the 1980s [20,21], the RCWA method has become a very mature algorithm after continuous development and optimization [22–27]. Its computing speed and numerical stability are constantly improved, and more and more application scenarios are also available. In addition, commercial software based on this method also continuously emerged, such as Ansys lumerical-RCWA [28] and MC grating [29]. However, there are not many plasmonic lithography models based on this method, especially the three-dimensional (3D) model building and how to improve the speed of 3D model calculation. Because with the development of plasmonic lithography technology, computational lithography technologies for plasmonic lithography such as source optimization (SO), mask optimization (MO), source mask optimization (SMO) and optical proximity correction (OPC), etc., will also be developed. These technologies will further improve the imaging level of plasmonic lithography and improve the resolution of plasmonic lithography, which is an important tool to promote the industrialization of the lithography technology. Figure 2 is a typical mask optimization example for a two-dimensional mask pattern [1]. However, computational lithography requires repeated invocation of the imaging model, as shown in Fig. 2(a), so the rapid calculation of the imaging model is highly required. Therefore, it is necessary to establish a fast and accurate optical imaging model of plasmonic lithography.

 figure: Fig. 2.

Fig. 2. A mask optimization schematic diagram: (a) the flow chart, (b) comparison of the effect before and after correction.

Download Full Size | PDF

This paper is divided into the following four parts: The first part describes the model parameters of the typical imaging structure of plasmonic lithography and the main principle and status of the existing RCWA algorithm. The second part is the improvement and innovation of the existing RCWA algorithm for the needs of computational lithography, that is, deducing the formula of calculating the light field inside the structure and innovatively proposing to integrate, store and invoke the scattering matrix to improve the calculation speed. The third part takes commercial software as the comparison object to carry on the comparison verification. The fourth part summarizes and looks forward.

2. Principles

2.1 Model parameters

Since the concept of surface plasmon was proposed and introduced into the field of optical lithography around the 2000s, plasmonic lithography has been continuously developed and has become an alternative optical lithography technology. Plasmonic lithography mainly includes plasmonic imaging lithography, plasmonic interference lithography and plasmonic direct writing lithography [2–5], in which imaging lithography and interference lithography have a mask structure, while direct writing lithography usually has mask-less structure. The imaging system of imaging lithography and interference lithography can be a superlens structure of a single layer of metal, as shown in Fig. 3(a); It can also be a hyperbolic multilayer metamaterial (HMM) structure with alternating metal/dielectric, as shown in Fig. 3(b) [30]. Imaging lithography is to achieve 1:1 imaging of mask pattern and photoresist image, and interference lithography is to use interference effect to achieve miniaturized imaging with 2:1 reduction, 4:1 reduction, … [30]. Direct writing lithography uses the local ability of the nanoantenna to the light field to strongly limit the excited local surface plasmons to the tip of the structure and form subwavelength light spots through the local surface plasmon resonance (SPR), thus achieving the purpose of pattern transfer. The schematic diagram is shown in Fig. 3(c) [31]. This paper focuses on the modeling of plasmonic lithography with mask structure, namely, plasmonic imaging lithography and plasmonic interference lithography.

 figure: Fig. 3.

Fig. 3. Schematic diagrams of three typical imaging structures for plasmonic lithography: (a) single-layer Ag film superlens imaging lithography structure, (b) hyperbolic multilayer film interference lithography structure, and (c) direct writing lithography system based on bowtie aperture. Reproduced from: (a–b), [30], Copyright © 2018, De Gruyter; (c), [31], Copyright © 2006, American Chemical Society.

Download Full Size | PDF

It can be seen from Fig. 1(a) and Fig. 3(a) and (b) that the entire optical system is stacked with metal or dielectric films of nanoscale thickness, and the arrangement also has a certain periodicity. The metal Cr is used as the shading part to form the mask layer. The mask pattern can be one-dimensional periodic pattern, such as the typical line/space pattern, or two-dimensional periodic pattern, such as regular holes and general irregular patterns. In general, the mid-plane of the photoresist PR layer is selected as the image plane of optical imaging. Considering the evanescent wave property, the thickness of the photoresist layer is also tens of nanometers. In addition, a metal film is often placed behind the photoresist layer to further improve the imaging effect of the photoresist layer through the metal reflection resonance effect [32].

In summary, the target modeling structure is a 3D multi-layer structure with several nanoscale thickness films arranged in the z direction of light wave propagation which is microscale in the horizontal x-y direction. The mask layer is the non-uniform layer with periodic pattern, and the other layers are uniform layer. The imaging plane is inside the structure rather than on the exit plane.

2.2 Introduction to the existing RCWA algorithm

The model established in this paper is mainly based on the 3D RCWA method, which is mostly used to calculate the reflectivity of the incident plane and the transmittance of the exit plane of the structure. The incident plane is the interface between the incident medium region and the first layer of the structure, and the exit plane is the interface between the transmitted medium region and the last layer of the structure. The incident region and the transmitted region are both semi-infinite spaces. In Ref. [19], the reflectance and transmittance of 3D structures, namely the electromagnetic field distribution of the incident plane and the exit plane, have been solved, and the scattering matrix has been further improved to improve the computational efficiency and numerical stability. The work of this paper is mainly based on the improvement and innovation of Ref. [19].

According to the formulas given in Ref. [19], the tangential electric field components ${s_{x,i}}$ and ${s_{y,i}}$, and magnetic field components ${u_{x,i}}$ and ${u_{y,i}}$ in Fourier space at different z coordinates of layer i of the structure can be expressed in matrix form as:

$${\Psi _i}({\tilde{z}} )= \left[ {\begin{array}{l} {{s_{x,i}}({\tilde{z}} )}\\ {{s_{y,i}}({\tilde{z}} )}\\ {{u_{x,i}}({\tilde{z}} )}\\ {{u_{y,i}}({\tilde{z}} )} \end{array}} \right] = \left[ {\begin{array}{cc} {{W_i}}&{{W_i}}\\ { - {V_i}}&{{V_i}} \end{array}} \right]\left[ {\begin{array}{cc} {{e^{ - {\lambda_i}\tilde{z}}}}&0\\ 0&{{e^{{\lambda_i}\tilde{z}}}} \end{array}} \right]\left[ {\begin{array}{c} {c_i^ + }\\ {c_i^ - } \end{array}} \right],$$
where light waves travel in the z direction. $\tilde{z} = {k_0}z$ and ${k_0} = \frac{{2\pi }}{\lambda }$ is the vacuum wave number. Both ${W_i}$ and ${V_i}$ are matrices that contain eigenvalue information of the i layer. The matrix ${W_i}$ describes the eigen-modes of the electric fields and ${V_i}$ describes the eigen-modes of the magnetic fields. ${\lambda _i}$ is the eigenvalue matrix of this layer and the exponential terms ${e^{ - {\lambda _i}\tilde{z}}}$ and ${e^{{\lambda _i}\tilde{z}}}$ describe forward and backward propagation through the layer respectively. $c_i^ + $ and $c_i^ - $ are the forward and back propagation coefficients in this layer, respectively. The detailed formula derivation can be seen in Ref. [19], which will not be repeated here. ${s_{x,i}}({\tilde{z}} )$, ${s_{y,i}}({\tilde{z}} )$, ${u_{x,i}}({\tilde{z}} )$ and ${u_{y,i}}({\tilde{z}} )$ respectively derived from the Fourier series expansion:
$$\begin{array}{l}{E_{x,i}}({x,y,z} )= \mathop \sum \limits_{m ={-} \infty }^\infty \mathop \sum \limits_{n ={-} \infty }^\infty {s_{x,i}}({m,n;z} )\cdot {e^{ - j[{{k_x}({m,n} )x + {k_y}({m,n} )y} ]}}\\{E_{y,i}}({x,y,z} )= \mathop \sum \limits_{m ={-} \infty }^\infty \mathop \sum \limits_{n ={-} \infty }^\infty {s_{y,i}}({m,n;z} )\cdot {e^{ - j[{{k_x}({m,n} )x + {k_y}({m,n} )y} ]}}\\{H_{x,i}}({x,y,z} )= \mathop \sum \limits_{m ={-} \infty }^\infty \mathop \sum \limits_{n ={-} \infty }^\infty {u_{x,i}}({m,n;z} )\cdot {e^{ - j[{{k_x}({m,n} )x + {k_y}({m,n} )y} ]}}\\{H_{y,i}}({x,y,z} )= \mathop \sum \limits_{m ={-} \infty }^\infty \mathop \sum \limits_{n ={-} \infty }^\infty {u_{y,i}}({m,n;z} )\cdot {e^{ - j[{{k_x}({m,n} )x + {k_y}({m,n} )y} ]}}\end{array}$$
where m = -∞, …, -1, 0, 1, …, ∞ and n = -∞, …, -1, 0, 1, …, ∞ are expansion orders of ${k_x}$ and ${k_y}$, respectively.

Equation (1) can be used to calculate the field throughout a layer that is uniform in the z direction. All real devices are composed of multiple layers and eigenmodes must be calculated separately for each layer. To connect the layers, boundary conditions are enforced by equating the tangential components of the fields on either side of the interfaces. In this manner, propagation through the entire device can be described rigorously. To solve the boundary condition problem efficiently, scattering matrices have arisen as the most popular method and solve the boundary condition problem one layer at a time instead of all layers simultaneously. The scattering matrix ${S^{(i )}}$ of layer i of the structure is defined as:

$$\left[ {\begin{array}{c} {c_{i - 1}^{^{\prime} - }}\\ {c_{i + 1}^{^{\prime} + }} \end{array}} \right] = {S^{(i )}}\left[ {\begin{array}{c} {c_{i - 1}^{^{\prime} + }}\\ {c_{i + 1}^{^{\prime} - }} \end{array}} \right],\; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; {S^{(i )}} = \left[ {\begin{array}{cc} {S_{11}^{(i )}}&{S_{12}^{(i )}}\\ {S_{21}^{(i )}}&{S_{22}^{(i )}} \end{array}} \right].$$

The schematic diagram is shown in Fig. 4(a). The column vectors $c_{i - 1}^{^{\prime} + }$, $c_{i - 1}^{^{\prime} - }$, $c_{i + 1}^{^{\prime} + }$ and $c_{i + 1}^{^{\prime} - }$ contain the mode coefficients of the field immediately outside of the ith layer in both forward and backward directions. Their subscripts indicate which side of the layer is being described while the signs in their superscripts indicate the direction of propagation.

 figure: Fig. 4.

Fig. 4. Schematic diagram of scattering matrix and propagation coefficient: (a) single scattering matrix of layer i, (b) global scattering matrix of the whole structure.

Download Full Size | PDF

By defining the spacing layer of zero thickness for wrapping and the Redheffer Star Product operation symbol ${\otimes} $, the scattering matrix ${S^{({global} )}}$ of the entire structure described in Fig. 4(b) connecting the incident substrate ref and the transmission substrate trn can be obtained, which is defined as:

$${S^{({global} )}} = {S^{({ref} )}} \otimes {S^{({mask} )}} \otimes {S^{({PMMA} )}} \otimes {S^{({Al} )}} \otimes {S^{({Si{O_2}} )}} \otimes \ldots \otimes {S^{({PR} )}} \otimes {S^{({trn} )}}.$$

The scattering matrix superscripts characterize the region respectively, for example ref and trn represent the incident substrate region and the transmission substrate region respectively. The scattering matrix ${S^{({global} )}}$ is expressed as:

$$\left[ {\begin{array}{c} {{c_{ref}}}\\ {{c_{trn}}} \end{array}} \right] = {S^{({global} )}}\left[ {\begin{array}{c} {{c_{inc}}}\\ 0 \end{array}} \right],{\kern 1cm} {S^{({global} )}} = \left[ {\begin{array}{cc} {S_{11}^{({global} )}}&{S_{12}^{({global} )}}\\ {S_{21}^{({global} )}}&{S_{22}^{({global} )}} \end{array}} \right].$$

The schematic diagram is shown in Fig. 4(b). Based on this global scattering matrix, the Fourier spatial frequency distribution of the tangential components of the reflection field ${\left[ {\begin{array}{cc} {s_x^{ref}}&{s_y^{ref}} \end{array}} \right]^T}$ and the transmitted field ${\left[ {\begin{array}{cc} {s_x^{trn}}&{s_y^{trn}} \end{array}} \right]^T}$ can be obtained according to the spatial spectrum of the incident electric field ${\left[ {\begin{array}{cc} {s_x^{inc}}&{s_y^{inc}} \end{array}} \right]^T}$:

$$\begin{array}{l}\left[ {\begin{array}{l} {s_x^{ref}}\\ {s_y^{ref}} \end{array}} \right] = {W_{ref}}{c_{ref}} = {W_{ref}}S_{11}^{({global} )}{c_{inc}} = {W_{ref}}S_{11}^{({global} )}W_{ref}^{ - 1}\left[ {\begin{array}{c} {s_x^{inc}}\\ {s_y^{inc}} \end{array}} \right]\\\left[ {\begin{array}{c} {s_x^{trn}}\\ {s_y^{trn}} \end{array}} \right] = {W_{trn}}{c_{trn}} = {W_{trn}}S_{21}^{({global} )}{c_{inc}} = {W_{trn}}S_{21}^{({global} )}W_{ref}^{ - 1}\left[ {\begin{array}{c} {s_x^{inc}}\\ {s_y^{inc}} \end{array}} \right]. \end{array}$$

The z components $s_z^{ref}$ and $s_z^{trn}$ of the reflected and transmitted electric fields can be solved by using the electric field divergence equation $\nabla \cdot \vec{E} = 0$, one of Maxwell's equations. The specific expression is:

$$\begin{array}{l}s_z^{ref} ={-} {\widetilde {{K_{z,ref}}}^{ - 1}}({\widetilde {{K_x}}s_x^{ref} + \widetilde {{K_y}}s_y^{ref}} )\\s_z^{trn} ={-} {\widetilde {{K_{z,trn}}}^{ - 1}}({\widetilde {{K_x}}s_x^{trn} + \widetilde {{K_y}}s_y^{trn}} ),\end{array}$$
where $\widetilde {{K_x}} = {k_x}/{k_0}$, $\widetilde {{K_y}} = {k_y}/{k_0}$ and $\widetilde {{K_z}} = {k_z}/{k_0}$. Then the electric field distribution in the spatial frequency domain can be obtained by Fourier transform. At this point, the main framework of algorithms for solving the distribution of reflection field and transmission field represented by Ref. [19] has been introduced. Next is the algorithm improvement we made to calculate the distribution of the inside field of the structure and improve the computational efficiency for computational lithography.

3. Algorithm improvement and innovation

3.1 Solving the inside field distribution of the structure

As shown by the black dotted line in Fig. 1(a) or the red dotted line in Fig. 4(b), the target of our solution is the electric field distribution of the middle section of the photoresist layer, which belongs to the inside of the structure rather than the exit plane. Therefore, the first improvement is to derive the formula for solving the electric field distribution in the middle section of the photoresist layer. The derivation process is as follows:

Equation (1) provides the expressions of tangential electric field and magnetic field components in Fourier space at different z coordinates of layer i of the structure, then the electromagnetic field distribution of the middle section of the photoresist layer can be expressed as:

$${\Psi _{PR}}({PR/2} )= \left[ {\begin{array}{c} {{s_{x,PR}}({PR/2} )}\\ {{s_{y,PR}}({PR/2} )}\\ {{u_{x,PR}}({PR/2} )}\\ {{u_{y,PR}}({PR/2} )} \end{array}} \right] = \left[ {\begin{array}{cc} {{W_{PR}}}&{{W_{PR}}}\\ { - {V_{PR}}}&{{V_{PR}}} \end{array}} \right]\left[ {\begin{array}{cc} {{e^{ - {\lambda_{PR}}{k_0}{L_{PR}}/2}}}&0\\ 0&{{e^{{\lambda_{PR}}{k_0}{L_{PR}}/2}}} \end{array}} \right]\left[ {\begin{array}{cc} {c_{PR}^ + }\\ {c_{PR}^ - } \end{array}} \right],$$
where $\left[ {\begin{array}{c} {c_{PR}^ + }\\ {c_{PR}^ - } \end{array}} \right]$ is the unknown quantity to be solved, and the others are known quantities. Combined with the boundary condition that the tangential components of the electric fields on both sides of the interface are continuous, the field distribution before and after the exit plane (that is, the interface of the photoresist layer and the transmission substrate, the blue dashed line in Fig. 4(b)) is equal, and the expression is:
$$\left[ {\begin{array}{cc} {{W_{PR}}}&{{W_{PR}}}\\ { - {V_{PR}}}&{{V_{PR}}} \end{array}} \right]\left[ {\begin{array}{cc} {{e^{ - {\lambda_{PR}}{k_0}{L_{PR}}}}}&0\\ 0&{{e^{{\lambda_{PR}}{k_0}{L_{PR}}}}} \end{array}} \right]\left[ {\begin{array}{c} {c_{PR}^ + }\\ {c_{PR}^ - } \end{array}} \right] = \left[ {\begin{array}{cc} {{W_{trn}}}&{{W_{trn}}}\\ { - {V_{trn}}}&{{V_{trn}}} \end{array}} \right]\left[ {\begin{array}{c} {{c_{trn}}}\\ 0 \end{array}} \right].$$

According to Eq. (9), the expression of $\left[ {\begin{array}{c} {c_{PR}^ + }\\ {c_{PR}^ - } \end{array}} \right]$ can be obtained, which can be substituted into Eq. (8) to obtain the expression of the tangential component of the electromagnetic field in the middle section of the photoresist layer:

$$\scalebox{0.7}{$\displaystyle\begin{array}{cc}{\Psi _{PR}}\left( {\frac{{PR}}{2}} \right) = \left[ {\begin{array}{cc} {{s_{x,PR}}\left( {\frac{{PR}}{2}} \right)}\\ {{s_{y,PR}}\left( {\frac{{PR}}{2}} \right)}\\ {{u_{x,PR}}\left( {\frac{{PR}}{2}} \right)}\\ {{u_{y,PR}}\left( {\frac{{PR}}{2}} \right)} \end{array}} \right]\\= \; \left[ {\begin{array}{cc} {{W_{PR}}}&{{W_{PR}}}\\ { - {V_{PR}}}&{{V_{PR}}} \end{array}} \right]{\left[ {\begin{array}{cc} {{e^{ - \frac{{{\lambda_{PR}}{k_0}{L_{PR}}}}{2}}}}&0\\ 0&{{e^{\frac{{{\lambda_{PR}}{k_0}{L_{PR}}}}{2}}}} \end{array}} \right]^{ - 1}}{\left[ {\begin{array}{cc} {{W_{PR}}}&{{W_{PR}}}\\ { - {V_{PR}}}&{{V_{PR}}} \end{array}} \right]^{ - 1}}\left[ {\begin{array}{cc} {{W_{trn}}}&{{W_{trn}}}\\ { - {V_{trn}}}&{{V_{trn}}} \end{array}} \right]\left[ {\begin{array}{c} {{c_{trn}}}\\ 0 \end{array}} \right]\\= \left[ {\begin{array}{cc} {{W_{PR}}}&{{W_{PR}}}\\ { - {V_{PR}}}&{{V_{PR}}} \end{array}} \right]{\left[ {\begin{array}{cc} {{e^{ - \frac{{{\lambda_{PR}}{k_0}{L_{PR}}}}{2}}}}&0\\ 0&{{e^{\frac{{{\lambda_{PR}}{k_0}{L_{PR}}}}{2}}}} \end{array}} \right]^{ - 1}}{\left[ {\begin{array}{cc} {{W_{PR}}}&{{W_{PR}}}\\ { - {V_{PR}}}&{{V_{PR}}} \end{array}} \right]^{ - 1}}\left[ {\begin{array}{cc} {{W_{trn}}}&{{W_{trn}}}\\ { - {V_{trn}}}&{{V_{trn}}} \end{array}} \right]\left[ {\begin{array}{c} {S_{21}^{({global} )}{c_{inc}}}\\ 0 \end{array}} \right]\\= \left[ {\begin{array}{cc} {{W_{PR}}}&{{W_{PR}}}\\ { - {V_{PR}}}&{{V_{PR}}} \end{array}} \right]{\left[ {\begin{array}{cc} {{e^{ - \frac{{{\lambda_{PR}}{k_0}{L_{PR}}}}{2}}}}&0\\ 0&{{e^{\frac{{{\lambda_{PR}}{k_0}{L_{PR}}}}{2}}}} \end{array}} \right]^{ - 1}}{\left[ {\begin{array}{cc} {{W_{PR}}}&{{W_{PR}}}\\ { - {V_{PR}}}&{{V_{PR}}} \end{array}} \right]^{ - 1}}\left[ {\begin{array}{cc} {{W_{trn}}}&{{W_{trn}}}\\ { - {V_{trn}}}&{{V_{trn}}} \end{array}} \right]\left[ {\begin{array}{c} {S_{21}^{({global} )}W_{ref}^{ - 1}\left[ {\begin{array}{c} {s_x^{inc}}\\ {s_y^{inc}} \end{array}} \right]}\\ 0 \end{array}} \right]\\= [{{A_{PR}}} ]\left[ {\begin{array}{*{20}{c}} {S_{21}^{({global} )}W_{ref}^{ - 1}\left[ {\begin{array}{*{20}{c}} {s_x^{inc}}\\ {s_y^{inc}} \end{array}} \right]}\\ 0 \end{array}} \right].\end{array}$}$$

The matrix $[{{A_{PR}}} ]= \left[ {\begin{smallmatrix} {{W_{PR}}}&{{W_{PR}}}\\ { - {V_{PR}}}&{{V_{PR}}} \end{smallmatrix}} \right]{\left[ {\begin{smallmatrix} {{e^{ - \frac{{{\lambda_{PR}}{k_0}{L_{PR}}}}{2}}}}&0\\ 0&{{e^{\frac{{{\lambda_{PR}}{k_0}{L_{PR}}}}{2}}}} \end{smallmatrix}} \right]^{ - 1}}{\left[ {\begin{smallmatrix} {{W_{PR}}}&{{W_{PR}}}\\ { - {V_{PR}}}&{{V_{PR}}} \end{smallmatrix}} \right]^{ - 1}}\left[ {\begin{smallmatrix} {{W_{trn}}}&{{W_{trn}}}\\ { - {V_{trn}}}&{{V_{trn}}} \end{smallmatrix}} \right]$ is characterized as a whole matrix, and is only determined by the geometrical parameters and material parameters of the photoresist layer. The superscript -1 represents the matrix inversion operation. Considering that there are both forward and back propagating fields in the structure, the z component ${s_{z,PR}}({PR/2} )$ of the electric field is solved by using another Maxwell equation, the magnetic field curl equation $\nabla \times \vec{H}{\varepsilon _r}\vec{E}$. The specific expression is 

$${s_{z,PR}}({PR/2} )= \frac{{\left( {\widetilde {{K_x}}{u_{y,PR}}\left( {\frac{{PR}}{2}} \right) - \widetilde {{K_y}}{u_{x,PR}}\left( {\frac{{PR}}{2}} \right)} \right)}}{{j{\varepsilon _{PR}}}},$$
where j is a purely imaginary number. Similarly, Fourier transform of the spatial frequency domain electric field distribution can obtain the spatial domain electric field distribution, then the light intensity distribution of the middle section of the photoresist layer can be expressed as:
$${I_{\frac{{PR}}{2}}}({x,y} )= {\left|{F{T^2}\left\{ {{s_{x,PR}}\left( {\frac{{PR}}{2}} \right)} \right\}} \right|^2} + {\left|{F{T^2}\left\{ {{s_{y,PR}}\left( {\frac{{PR}}{2}} \right)} \right\}} \right|^2} + {\left|{F{T^2}\left\{ {{s_{z,PR}}\left( {\frac{{PR}}{2}} \right)} \right\}} \right|^2},$$
where $F{T^2}\{{\; \cdot \; } \}$ represents the two-dimensional Fourier transform and $|{\; \cdot \; } |$. represents the modulo operation.

3.2 Engineering improvement for computational lithography requirements – integration, storage and invocation of scattering matrices

The second innovation we have made is closely focused on the theme of computational lithography, with the aim of increasing the speed of photoresist aerial image calculation. From the perspective of practical engineering applications, the actual practice of mask optimization often changes the mask pattern shape, while the lighting mode and imaging structure (the multilayer film between the mask and the photoresist) are relatively fixed. For optical imaging from mask pattern to photoresist, the main calculation time of RCWA algorithm should be spent on calculating the scattering matrix of the mask layer when calculating different patterns, because the change of the mask pattern will cause the change of the scattering matrix of the mask layer. As an imaging system, the middle multilayer structure generally does not change, so the scattering matrix of the middle multilayer film does not change. Therefore, it can be stored and called after the first calculation, so as to save time. Specifically, the Eq. (4) is rewritten as :

$$\begin{aligned}{S^{({global} )}}({MASK} )&= {S^{({ref} )}} \otimes {S^{({mask} )}} \otimes {S^{({PMMA} )}} \otimes {S^{({Al} )}} \otimes {S^{({Si{O_2}} )}} \otimes \ldots \otimes {S^{({PR} )}} \otimes {S^{({trn} )}}\\&= {S^{({ref} )}} \otimes {S^{({mask} )}} \otimes {S^{({TRN} )}},\end{aligned}$$
where ${S^{({TRN} )}} = {S^{({PMMA} )}} \otimes {S^{({Al} )}} \otimes {S^{({Si{O_2}} )}} \otimes \ldots \otimes {S^{({PR} )}} \otimes {S^{({trn} )}}$ is regarded as a whole scattering matrix. And several scattering matrices contained in it need only be solved once and calculated in this way to obtain ${S^{({TRN} )}}$, then stored and called. Similarly, ${S^{({ref} )}}$ also only needs to be calculated once and then stored and called. When mask optimization is performed, the mask pattern changes and the corresponding ${S^{({mask} )}}$ changes with it, but ${S^{({ref} )}}$ and ${S^{({TRN} )}}$ do not change. Therefore, the main computational effort should be spent on the calculation of ${S^{({mask} )}}$ when calculating ${S^{({global} )}}({MASK} )$. At this time, when the mask pattern changes, the expression of the tangential component of the electromagnetic field in the middle section of the photoresist layer can be expressed as :
$$\begin{aligned}{\Psi _{PR}}\left( {\frac{{PR}}{2},MASK} \right) &= \left[ {\begin{array}{c} {{s_{x,PR}}\left( {\frac{{PR}}{2},MASK} \right)}\\ {{s_{y,PR}}\left( {\frac{{PR}}{2},MASK} \right)}\\ {{u_{x,PR}}\left( {\frac{{PR}}{2},MASK} \right)}\\ {{u_{y,PR}}\left( {\frac{{PR}}{2},MASK} \right)} \end{array}} \right] = [{{A_{PR}}} ]\left[ {\begin{array}{c} {S_{21}^{({global} )}({MASK} )W_{ref}^{ - 1}\left[ {\begin{array}{c} {s_x^{inc}}\\ {s_y^{inc}} \end{array}} \right]}\\ 0 \end{array}} \right]\\&= [{{A_{PR}}} ]\left[ {\begin{array}{c} {{{[{{S^{({ref} )}} \otimes {S^{({mask} )}} \otimes {S^{({TRN} )}}} ]}_{21}}W_{ref}^{ - 1}\left[ {\begin{array}{c} {s_x^{inc}}\\ {s_y^{inc}} \end{array}} \right]}\\ 0 \end{array}} \right].\end{aligned}$$

The solution of z component of electric field and light intensity can be obtained by operating Eqs. (11), (12). In this way, it is only necessary to calculate ${S^{({mask} )}}$ and then perform a few simple matrix operations to get the light intensity distribution of the new pattern. This will significantly reduce the calculation time, especially the more layers of the structure behind the mask, the more time saved.

In the following, we will show the results of calculating the aerial image intensity distribution of the photoresist middle section based on the improved RCWA algorithm, and compare them with some other numerical software to show the advantages of this improved algorithm in modeling.

4. Comparison verification

It is reported that the latest version R1.3 of ANSYS Lumerical has software with GUI based on the RCWA algorithm, and its accuracy has been tested [33]. Therefore, we use this commercial software as a comparison object. Through the test of several typical patterns, the improved RCWA algorithm is verified to improve the calculation speed compared with the software. At the same time, the error of the two calculation results will also be presented to show the accuracy of our improved algorithm. Our improved algorithm and the commercial software to be compared are both running on the same server (Intel Xeon Gold 6226R CPU @ 2,90 GHz, with 256 GB RAM).

The superlens-based imaging lithography device shown in Fig. 1(a) is employed. The wavelength λ is 365 nm, and the substrate of the mask is quartz glass, with permittivity of 1.46. Along the positive direction of light propagation, the mask Cr grating layer, PMMA spacer layer, Ag film and photoresist PR are sequentially distributed. The corresponding thickness of each layer is 50 nm, 50 nm, 20 nm and 40 nm. And the corresponding permittivities of each layer is -8.55 + 8.96i, 2.25, -2.4 + 0.25i and 2.56, respectively [30]. The transmission substrate is the reflective Ag layer with permittivity of -2.4 + 0.25i. Typical patterns selected are square holes and general T-shaped general patterns.

4.1 Square holes

Taking into account the resolution level of the superlens at 365 nm wavelength, the square hole pattern period is set to 400 nm and the through-light hole CD is set to 200 nm, as shown in Fig. 5. The illumination mode is X-polarized light with normal incidence, and the imaging plane was the middle plane of the photoresist layer (Fig. 1(a) black dashed line). The sampling point is 1 nm × 1 nm, and the expansion order is ${k_x}$ = ${k_y}$ = [-15, -14, …, 0, …, 14, 15]. The convergence analysis of the expanded order is shown in section 4.3.

 figure: Fig. 5.

Fig. 5. Schematic diagram of the mask pattern (square hole), where blue is transparent (PMMA) and yellow is opaque (Cr).

Download Full Size | PDF

Figure 6 shows the calculation results of the improved RCWA algorithm and ANSYS software respectively. It has been explained in Section 3.1 that Ex, Ey and Ez need to be calculated separately. Therefore, Fig. 6 (a)–(c) are normalized amplitude components Ex, Ey and Ez calculated by the improved RCWA algorithm, while Fig. 6 (d)–(f) are normalized amplitude components Ex, Ey and Ez calculated by ANSYS software, respectively. From a qualitative point of view, the results obtained by the two methods are very close.

 figure: Fig. 6.

Fig. 6. Comparison of the calculated results of the improved RCWA algorithm and ANSYS software, both of which are normalized amplitude components:(a) Ex of improved RCWA algorithm, (b) Ey of improved RCWA algorithm, (c) Ez of improved RCWA algorithm;(d) Ex of ANSYS software, (e) Ey of ANSYS software, and (f) Ez of ANSYS software.

Download Full Size | PDF

Furthermore, according to the amplitude components, the intensity distributions calculated by the two methods are obtained respectively. Figure 7(a) shows the normalized light intensity distribution of the improved RCWA algorithm, and Fig. 7(b) shows the normalized light intensity distribution of ANSYS software. The distribution error of light intensity obtained by the two methods is shown in Fig. 7 (c), and the maximum error is about 0.08 (V/m)2, and the root mean square error (RMSE) is 0.0337 (V/m)2. The error and RMSE are calculated in Eq. (15):

$$RMSE = \sqrt {\mathop \sum \limits_{x,y} \frac{{error{{({x,y} )}^2}}}{{{N_x} \times {N_y}}}} ,\textrm{e}rror({x,y} )= \,{I_{M - RCWA}}({x,y} )- {I_{ANSYS}}({x,y} ). $$

In addition, at the expansion order of ${k_x}$ = ${k_y}$= [-15, -14, …, 0, …, 14, 15], the calculation time required by ANSYS software is about 165 seconds, while the first calculation of our improved RCWA algorithm takes only 53 seconds, and the speed is increased by 2.1 times. This might be mainly owing to the difference in the way the two calculate the scattering matrix of each layer and the difference in algorithm adaptation. Our improved RCWA algorithm is built in MATLAB, which is more suitable for matrix operations, so it is slightly faster than ANSYS software without the separation and integration of scattering matrices. In addition, when the method of Eq. (14) is used to integrate, store and call other matrices required for calculation other than ${S^{({mask} )}}$, only change ${S^{({mask} )}}$ and make it participate in the calculation, the operation time is reduced to 20 seconds. At this point, the computing speed is increased by 7.25 times.

 figure: Fig. 7.

Fig. 7. Comparison between the normalized light intensity calculation results of the improved RCWA algorithm and ANSYS software: (a) the normalized light intensity of the improved RCWA algorithm, (b) the normalized light intensity of the ANSYS software, and (c) the error distribution of the two.

Download Full Size | PDF

4.2 General T-shaped pattern

The imaging device remains unchanged, and the general T-shaped pattern period is increased and set to 800 nm, while the light transmission area is a T-shaped pattern composed of rectangular holes of 600 nm × 200 nm, as shown in Fig. 8. The illumination mode is still the normal incident X-polarized light, and the imaging plane is still the middle plane of the photoresist layer (Fig. 1(a) black dashed line). The sampling point is still a 1 nm × 1 nm square, and the expansion order is ${k_x}$ = ${k_y}$= [-20, -19, …, 0, …, 19, 20]. Similarly, the hierarchical convergence analysis is presented in section 4.3.

 figure: Fig. 8.

Fig. 8. Schematic#diagram of the mask pattern (general T-shaped pattern), where blue is transparent (PMMA) and yellow is opaque (Cr).

Download Full Size | PDF

Similar to Section 4.1, Fig. 9 and Fig. 10 show the amplitude, light intensity and light intensity error distributions of each component of the electric field calculated by the two methods respectively. From the qualitative point of view, the results of the two methods are still very close. And from the quantitative perspective, the maximum error of light intensity distribution is about 0.12 (V/m)2, and the root mean square error (RMSE) is 0.0385 (V/m)2. In addition, at the expansion order of ${k_x}$ = ${k_y}$= [-20, -19, …, 0, …, 19, 20], the calculation time required by ANSYS software is about 780 seconds, while the first calculation of our improved RCWA algorithm only takes 222 seconds, and the speed is increased by 2.5 times. In addition, when the method of Eq. (14) in Section 3.2 is adopted, that is, other matrices required for calculation other than ${S^{({mask} )}}$ are integrated, stored and then invoked, and only ${S^{({mask} )}}$ is changed and allowed to participate in the calculation, then the operation time is reduced to 83 seconds. At this point, the computing speed is increased by 8.2 times.

 figure: Fig. 9.

Fig. 9. Comparison of the calculated results of the improved RCWA algorithm and ANSYS software, both of which are normalized amplitude components:(a) Ex of improved RCWA algorithm, (b) Ey of improved RCWA algorithm, (c) Ez of improved RCWA algorithm;(d) Ex of ANSYS software, (e) Ey of ANSYS software, and (f) Ez of ANSYS software.

Download Full Size | PDF

 figure: Fig. 10.

Fig. 10. Comparison between the normalized light intensity calculation results of the improved RCWA algorithm and ANSYS software: (a) the normalized light intensity of the improved RCWA algorithm, (b) the normalized light intensity of the ANSYS software, and (c) the error distribution of the two.

Download Full Size | PDF

4.3 Convergence analysis of the expanded order

In the evaluation of RCWA algorithm, the relevant parameters include running time and accuracy, and the most relevant physical quantity is the expansion order. Of course, the number of layers of the structure will also affect the calculation time, and the more the number of layers, the longer the calculation time. In the whole verification process, the same imaging structure is adopted, so the number of layers of the structure is also the same, so as to ensure that the comparison of results will not be unfair due to the different influencing factors of the number of layers of the structure.

At normal incidence, the expansion order can be expressed as

$$\begin{array}{l}{k_x} = m\frac{{2\pi }}{{{D_x}}}\;\;\;m = - \infty , \ldots \; , - 2, - 1,\; 0,\; 1,\; 2, \ldots \; \; ,\; \infty,\\{k_y} = n\frac{{2\pi }}{{{D_y}}},\;\;\;n = - \infty , \ldots \; , - 2, - 1,\; 0,\; 1,\; 2, \ldots \; \; ,\; \infty,\end{array}$$
where m and n are the number of expansion orders, respectively, and ${D_x}$ and ${D_y}$ are the periods in the x and y directions, respectively. Generally speaking, the more expansion orders, the higher the accuracy of the calculation, but at the same time, the more time it takes to calculate. Therefore, from the perspective of practical application, a limited number of expansion orders are often used, so it is also called truncation orders. In this way, the calculation accuracy can be guaranteed and the calculation speed can be fast, so it is necessary to choose the appropriate expansion order.

In plasmonic lithography, since the imaging range is localized in the near field and evanescent wave is also involved in imaging through the excitation of SPP, the frequency band of the expansion order should contain both transmission wave and evanescent wave. References [15,30,32] has calculated the optical transfer function (OTF) curve of superlens structure, as shown in Fig. 11. The curve corresponding to 0 to $n{k_0}$ is the transmission wave band, and the range greater than $n{k_0}$ is the evanescent wave band. The formant after $n{k_0}$ is the spatial frequency ${k_{spp}}$ corresponding to the SPP. After exceeding ${k_{spp}}$, the transmission efficiency decreases with the increase of spatial frequency, so that it can be basically ignored at $10{k_0}$. Moreover, the diffraction order of the mask also decreases with the increase of the spatial frequency. Therefore, it is appropriate to choose the truncation order near $10{k_0}$. At this time, the calculation accuracy and calculation speed can be balanced.

 figure: Fig. 11.

Fig. 11. The OTF curve of the superlens structure.

Download Full Size | PDF

According to Eq. (16), in order to achieve the same frequency range, the larger the period of the pattern, the more expansion orders are required. The pattern period of section 4.1 is 400 nm and the maximum order is 15, so the maximum spatial frequency that can be achieved is about ${k_{x,max}} = {k_{y,max}} = m\frac{{2\pi }}{{{D_x}}} = n\frac{{2\pi }}{{{D_y}}} = 15\frac{{2\pi }}{{400nm}} = 15\frac{{2\pi /400nm}}{{2\pi /365nm}}{k_0} = 13.69{k_0}$; The pattern period of section 4.2 is 800 nm and the maximum order is 20, so the maximum spatial frequency that can be achieved is about ${k_{x,max}} = {k_{y,max}} = m\frac{{2\pi }}{{{D_x}}} = n\frac{{2\pi }}{{{D_y}}} = 20\frac{{2\pi }}{{800nm}} = 20\frac{{2\pi /800nm}}{{2\pi /365nm}}{k_0} = 9.13{k_0}$. In order to show that different expansion orders affect the calculation time, different expansion orders are used in sections 4.1 and 4.2 under the condition that convergence is guaranteed. The expansion order of section 4.2 is higher than that of section 4.1, so the calculation time is also longer. The results of comparison and analysis show the accuracy of the calculation results with the above expansion orders.

In addition, in order to fully verify the universality of the proposed model, one-dimensional periodic patterns and two-dimensional patterns under partially coherent illumination conditions are also compared and verified respectively. For specific data, see the Supplement 1.

5. Conclusion

In this paper, the principle and characteristics of plasmonic lithography are reminded first, and the importance of establishing plasmonic lithography imaging model and developing computational lithography is pointed out. The implement of these works is inseparable from the calculation, especially the computational lithography needs to repeatedly call the imaging model for calculation. Therefore, the establishment of a model based on fast and high-precision algorithms is very urgent.

According to the characteristics of plasmonic lithography, the RCWA algorithm is a very suitable alternative algorithm. Next, the model parameters of typical imaging structures in plasmonic lithography and the principle and status of existing RCWA algorithm are described. Then, we focus on the improvement and innovation of the existing RCWA algorithm for the needs of computational lithography, that is, the formula for calculating the inside light field of the structure is derived and the innovation is proposed to integrate, store and call the scattering matrix to improve the calculation speed. Finally, the results are compared with commercial software for the two typical patterns. The results show that the two calculation results are very close, with RMSE less than 0.04 (V/m)2. The results obtained by our improved RCWA algorithm are slightly different from those obtained by ANSYS software. In addition to the inherent numerical error, the main reason is that the two methods of calculating the scattering matrix of each layer are different. RMSE is related to the specific value of light intensity, and the RMSE of the two is about two orders of magnitude smaller than the calculated specific result of light intensity, which is considered to be within the acceptable range, and the accuracy will be further improved in the future. Moreover, the calculation speed can be increased by more than 2 times during the first operation, and can be increased to about 8 times by integrating, storing and calling with the improved RCWA algorithm.

In summary, the improved RCWA algorithm is well suited to the needs of plasmonic lithography for computational lithography. The principle of the algorithm is based on the basic physical derivation and has passed the comparison and verification of commercial software. At the same time, the calculation speed has been greatly improved, which creates the conditions for the development of plasmonic computational lithography.

Funding

Special Project for Research and Development in Key areas of Guangdong Province (2022B0701180001); Guangzhou City Research and Development Program in Key Fields (202103020001); The construction of new research and development institutions (2019B090904015); A high-level innovation research institute from Guangdong Greater Bay Area Institute of Integrated Circuit and System (2019B090909006); Guangdong Province Research and Development Program in Key Fields (2021B0101280002); Fundamental Research Funds for the Central Universities (E2ET3801); University of Chinese Academy of Sciences Education Foundation[China] (118900M032).

Disclosures

The authors declare no conflicts of interest.

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.

Supplemental document

See Supplement 1 for supporting content.

References

1. X. Ma and G. R. Arce, Computational lithography (Wiley, 2010), pp. 5–7.

2. C. Wang, Z. Zhao, P. Gao, et al., “Surface plasmon lithography beyond the diffraction limit (in Chinese),” Chin. Sci. Bull. 61(6), 585–589 (2016). [CrossRef]  

3. C. Wang, W. Zhang, Z. Zhao, et al., “Plasmonic structures, materials and lenses for optical lithography beyond the diffraction limit: A review,” Micromachines 7(7), 118 (2016). [CrossRef]  

4. Y. Luo, X. Jiang, L. Liu, et al., “Recent Advances in Plasmonic Nanolithography,” Nano Lett. 10(1), 1–12 (2018). [CrossRef]  

5. F. Hong and R. Blaikie, “Plasmonic Lithography: Recent Progress,” Adv. Opt. Mater. 7(14), 1801653 (2019). [CrossRef]  

6. T. Ito, M. Ogino, T. Yamada, et al., “Fabrication of sub-100 nm Patterns using Near-field,” J. Photopolym. Sci. Technol. 18(3), 435–441 (2005). [CrossRef]  

7. P. Gao, N. Yao, C. Wang, et al., “Enhancing aspect profile of half-pitch 32 nm and 22 nm lithography with plasmonic cavity lens,” Appl. Phys. Lett. 106(9), 093110 (2015). [CrossRef]  

8. K. S. Yee, “Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. 14(3), 302–307 (1966). [CrossRef]  

9. J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. 114(2), 185–200 (1994). [CrossRef]  

10. Y. Zhu, “Multigrid Finite Element Methods for Electromagnetic Field Modeling,” IEEE Xplore, (2006).

11. N. Yadav and N. Kumar, “Analytical technique for electromagnetic field using finite element method,” International Journal of Engineering Trends & Technology 4(7), 3392 (2013).

12. A. A. Adetoyinbo and O. O. Adewole, “An algorithm for solving electromagnetic field equations by finite element method,” Research Journal of Applied Sciences 11, 1182 (2007).

13. C. P. Moore, R. J. Blaikie, and M. D. Arnold, “An improved transfer-matrix model for optical superlenses,” Opt. Express 17(16), 14260–14269 (2009). [CrossRef]  

14. Z. Zhao, Y. Luo, N. Yao, et al., “Modeling and experimental study of plasmonic lens imaging with resolution enhanced methods,” Opt. Express 24(24), 27115 (2016). [CrossRef]  

15. H. Ding, L. Liu, L. Dong, et al., “Study on forbidden pitch in plasmonic lithography: taking 365 nm wavelength thin silver film-based superlens imaging lithography as an example,” Opt. Express 30(19), 33869–33885 (2022). [CrossRef]  

16. H. Ding, L. Liu, Z. Li, et al., “Plasmonic lithography fast imaging model based on the decomposition machine learning method,” Opt. Express 31(1), 192–209 (2023). [CrossRef]  

17. N. Chateau and J. P. Hugonin, “Algorithm for the rigorous coupled-wave analysis of grating diffraction,” J. Opt. Soc. Am. A 11(4), 1321–1331 (1994). [CrossRef]  

18. D. M. Whittaker and I. S. Culshaw, “Scattering-matrix treatment of patterned multilayer photonic structures,” Phys. Rev. B 60(4), 2610–2618 (1999). [CrossRef]  

19. R. C. Rumpf, “Improved formulation of scattering matrices for semi-analytical methods that is consistent with convention,” PIER B 35, 241–261 (2011). [CrossRef]  

20. M. G. Moharam and T. K. Gaylord, “Diffraction analysis of dielectric surface-relief gratings,” J. Opt. Soc. Am. 72(10), 1385–1392 (1982). [CrossRef]  

21. M. G. Moharam, E. B. Grann, D. A. Pommet, et al., “Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings,” J. Opt. Soc. Am. A 12(5), 1068–1076 (1995). [CrossRef]  

22. S. Peng and G. M. Morris, “Efficient implementation of rigorous coupled-wave analysis for surface-relief gratings,” J. Opt. Soc. Am. A 12(5), 1087–1096 (1995). [CrossRef]  

23. P. Lalanne and G. M. Morris, “Highly improved convergence of the coupled-wave method for TM polarization,” J. Opt. Soc. Am. A 13(4), 779–784 (1996). [CrossRef]  

24. L. Li, “Use of Fourier series in the analysis of discontinuous periodic structures,” J. Opt. Soc. Am. A 13(9), 1870–1876 (1996). [CrossRef]  

25. N. Lyndin, O. Parriaux, and A.V. Tishchenko, “Modal analysis and suppression of the FMM instabilities in highly conductive gratings,” J. Opt. Soc. Am. A 24(12), 3781–3788 (2007). [CrossRef]  

26. V. Liu and S. Fan, “S4: A free electromagnetic solver for layered periodic structures,” Comput. Phys. Commun. 183(10), 2233–2244 (2012). [CrossRef]  

27. S. Spiridonov and A. A. Shcherbakov, “Reformulated Fourier Modal Method with improved near field computations,” J. Comput. Sci. 67, 101936 (2023). [CrossRef]  

28. Ansys Lumerical RCWA Simulation of Multilayer Periodic Photonic Structures (2023), https://www.ansys.com/products/photonics/ansys-lumerical-rcwa.

29. Modal and C Methods Grating Design and Analysis Software (2020) https://mcgrating.com/..

30. G. Liang, X. Chen, Q. Zhao, et al., “Achieving pattern uniformity in plasmonic lithography by spatial frequency selection,” Nanophotonics 7(1), 277–286 (2018). [CrossRef]  

31. L. Wang, S. M. Uppuluri, E. X. Jin, et al., “Nanolithography using high transmission nanoscale bowtie apertures,” Nano Lett. 6(3), 361–364 (2006). [CrossRef]  

32. T. Xu, L. Fang, J. Ma, et al., “Localizing surface plasmons with a metal-cladding superlens for projecting deep-subwavelength patterns,” Appl. Phys. B 97(1), 175–179 (2009). [CrossRef]  

33. Lumerical-Sub-wavelength-Model-Introduction-and-Data-Generation (2023) https://optics.ansys.com/hc/en-us/articles/8597760630163-Lumerical-Sub-wavelength-Model-Introduction-and-Data-Generation.

Supplementary Material (1)

NameDescription
Supplement 1       Supplement 1

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.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (11)
Fig. 1.
Fig. 1. Schematic diagram of imaging principle of plasmonic lithography and conventional projection optical lithography: (a) plasmonic lithography, (b) DUV lithography.

View in Article | Download Full Size | PDF

Fig. 2.
Fig. 2. A mask optimization schematic diagram: (a) the flow chart, (b) comparison of the effect before and after correction.

View in Article | Download Full Size | PDF

Fig. 3.
Fig. 3. Schematic diagrams of three typical imaging structures for plasmonic lithography: (a) single-layer Ag film superlens imaging lithography structure, (b) hyperbolic multilayer film interference lithography structure, and (c) direct writing lithography system based on bowtie aperture. Reproduced from: (a–b), [30], Copyright © 2018, De Gruyter; (c), [31], Copyright © 2006, American Chemical Society.

View in Article | Download Full Size | PDF

Fig. 4.
Fig. 4. Schematic diagram of scattering matrix and propagation coefficient: (a) single scattering matrix of layer i, (b) global scattering matrix of the whole structure.

View in Article | Download Full Size | PDF

Fig. 5.
Fig. 5. Schematic diagram of the mask pattern (square hole), where blue is transparent (PMMA) and yellow is opaque (Cr).

View in Article | Download Full Size | PDF

Fig. 6.
Fig. 6. Comparison of the calculated results of the improved RCWA algorithm and ANSYS software, both of which are normalized amplitude components:(a) Ex of improved RCWA algorithm, (b) Ey of improved RCWA algorithm, (c) Ez of improved RCWA algorithm;(d) Ex of ANSYS software, (e) Ey of ANSYS software, and (f) Ez of ANSYS software.

View in Article | Download Full Size | PDF

Fig. 7.
Fig. 7. Comparison between the normalized light intensity calculation results of the improved RCWA algorithm and ANSYS software: (a) the normalized light intensity of the improved RCWA algorithm, (b) the normalized light intensity of the ANSYS software, and (c) the error distribution of the two.

View in Article | Download Full Size | PDF

Fig. 8.
Fig. 8. Schematic#diagram of the mask pattern (general T-shaped pattern), where blue is transparent (PMMA) and yellow is opaque (Cr).

View in Article | Download Full Size | PDF

Fig. 9.
Fig. 9. Comparison of the calculated results of the improved RCWA algorithm and ANSYS software, both of which are normalized amplitude components:(a) Ex of improved RCWA algorithm, (b) Ey of improved RCWA algorithm, (c) Ez of improved RCWA algorithm;(d) Ex of ANSYS software, (e) Ey of ANSYS software, and (f) Ez of ANSYS software.

View in Article | Download Full Size | PDF

Fig. 10.
Fig. 10. Comparison between the normalized light intensity calculation results of the improved RCWA algorithm and ANSYS software: (a) the normalized light intensity of the improved RCWA algorithm, (b) the normalized light intensity of the ANSYS software, and (c) the error distribution of the two.

View in Article | Download Full Size | PDF

Fig. 11.
Fig. 11. The OTF curve of the superlens structure.

View in Article | Download Full Size | PDF

Equations (16)

Equations on this page are rendered with MathJax. Learn more.

$${\Psi _i}({\tilde{z}} )= \left[ {\begin{array}{l} {{s_{x,i}}({\tilde{z}} )}\\ {{s_{y,i}}({\tilde{z}} )}\\ {{u_{x,i}}({\tilde{z}} )}\\ {{u_{y,i}}({\tilde{z}} )} \end{array}} \right] = \left[ {\begin{array}{cc} {{W_i}}&{{W_i}}\\ { - {V_i}}&{{V_i}} \end{array}} \right]\left[ {\begin{array}{cc} {{e^{ - {\lambda_i}\tilde{z}}}}&0\\ 0&{{e^{{\lambda_i}\tilde{z}}}} \end{array}} \right]\left[ {\begin{array}{c} {c_i^ + }\\ {c_i^ - } \end{array}} \right],$$
$$\begin{array}{l}{E_{x,i}}({x,y,z} )= \mathop \sum \limits_{m ={-} \infty }^\infty \mathop \sum \limits_{n ={-} \infty }^\infty {s_{x,i}}({m,n;z} )\cdot {e^{ - j[{{k_x}({m,n} )x + {k_y}({m,n} )y} ]}}\\{E_{y,i}}({x,y,z} )= \mathop \sum \limits_{m ={-} \infty }^\infty \mathop \sum \limits_{n ={-} \infty }^\infty {s_{y,i}}({m,n;z} )\cdot {e^{ - j[{{k_x}({m,n} )x + {k_y}({m,n} )y} ]}}\\{H_{x,i}}({x,y,z} )= \mathop \sum \limits_{m ={-} \infty }^\infty \mathop \sum \limits_{n ={-} \infty }^\infty {u_{x,i}}({m,n;z} )\cdot {e^{ - j[{{k_x}({m,n} )x + {k_y}({m,n} )y} ]}}\\{H_{y,i}}({x,y,z} )= \mathop \sum \limits_{m ={-} \infty }^\infty \mathop \sum \limits_{n ={-} \infty }^\infty {u_{y,i}}({m,n;z} )\cdot {e^{ - j[{{k_x}({m,n} )x + {k_y}({m,n} )y} ]}}\end{array}$$
$$\left[ {\begin{array}{c} {c_{i - 1}^{^{\prime} - }}\\ {c_{i + 1}^{^{\prime} + }} \end{array}} \right] = {S^{(i )}}\left[ {\begin{array}{c} {c_{i - 1}^{^{\prime} + }}\\ {c_{i + 1}^{^{\prime} - }} \end{array}} \right],\; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; {S^{(i )}} = \left[ {\begin{array}{cc} {S_{11}^{(i )}}&{S_{12}^{(i )}}\\ {S_{21}^{(i )}}&{S_{22}^{(i )}} \end{array}} \right].$$
$${S^{({global} )}} = {S^{({ref} )}} \otimes {S^{({mask} )}} \otimes {S^{({PMMA} )}} \otimes {S^{({Al} )}} \otimes {S^{({Si{O_2}} )}} \otimes \ldots \otimes {S^{({PR} )}} \otimes {S^{({trn} )}}.$$
$$\left[ {\begin{array}{c} {{c_{ref}}}\\ {{c_{trn}}} \end{array}} \right] = {S^{({global} )}}\left[ {\begin{array}{c} {{c_{inc}}}\\ 0 \end{array}} \right],{\kern 1cm} {S^{({global} )}} = \left[ {\begin{array}{cc} {S_{11}^{({global} )}}&{S_{12}^{({global} )}}\\ {S_{21}^{({global} )}}&{S_{22}^{({global} )}} \end{array}} \right].$$
$$\begin{array}{l}\left[ {\begin{array}{l} {s_x^{ref}}\\ {s_y^{ref}} \end{array}} \right] = {W_{ref}}{c_{ref}} = {W_{ref}}S_{11}^{({global} )}{c_{inc}} = {W_{ref}}S_{11}^{({global} )}W_{ref}^{ - 1}\left[ {\begin{array}{c} {s_x^{inc}}\\ {s_y^{inc}} \end{array}} \right]\\\left[ {\begin{array}{c} {s_x^{trn}}\\ {s_y^{trn}} \end{array}} \right] = {W_{trn}}{c_{trn}} = {W_{trn}}S_{21}^{({global} )}{c_{inc}} = {W_{trn}}S_{21}^{({global} )}W_{ref}^{ - 1}\left[ {\begin{array}{c} {s_x^{inc}}\\ {s_y^{inc}} \end{array}} \right]. \end{array}$$
$$\begin{array}{l}s_z^{ref} ={-} {\widetilde {{K_{z,ref}}}^{ - 1}}({\widetilde {{K_x}}s_x^{ref} + \widetilde {{K_y}}s_y^{ref}} )\\s_z^{trn} ={-} {\widetilde {{K_{z,trn}}}^{ - 1}}({\widetilde {{K_x}}s_x^{trn} + \widetilde {{K_y}}s_y^{trn}} ),\end{array}$$
$${\Psi _{PR}}({PR/2} )= \left[ {\begin{array}{c} {{s_{x,PR}}({PR/2} )}\\ {{s_{y,PR}}({PR/2} )}\\ {{u_{x,PR}}({PR/2} )}\\ {{u_{y,PR}}({PR/2} )} \end{array}} \right] = \left[ {\begin{array}{cc} {{W_{PR}}}&{{W_{PR}}}\\ { - {V_{PR}}}&{{V_{PR}}} \end{array}} \right]\left[ {\begin{array}{cc} {{e^{ - {\lambda_{PR}}{k_0}{L_{PR}}/2}}}&0\\ 0&{{e^{{\lambda_{PR}}{k_0}{L_{PR}}/2}}} \end{array}} \right]\left[ {\begin{array}{cc} {c_{PR}^ + }\\ {c_{PR}^ - } \end{array}} \right],$$
$$\left[ {\begin{array}{cc} {{W_{PR}}}&{{W_{PR}}}\\ { - {V_{PR}}}&{{V_{PR}}} \end{array}} \right]\left[ {\begin{array}{cc} {{e^{ - {\lambda_{PR}}{k_0}{L_{PR}}}}}&0\\ 0&{{e^{{\lambda_{PR}}{k_0}{L_{PR}}}}} \end{array}} \right]\left[ {\begin{array}{c} {c_{PR}^ + }\\ {c_{PR}^ - } \end{array}} \right] = \left[ {\begin{array}{cc} {{W_{trn}}}&{{W_{trn}}}\\ { - {V_{trn}}}&{{V_{trn}}} \end{array}} \right]\left[ {\begin{array}{c} {{c_{trn}}}\\ 0 \end{array}} \right].$$
$$\scalebox{0.7}{$\displaystyle\begin{array}{cc}{\Psi _{PR}}\left( {\frac{{PR}}{2}} \right) = \left[ {\begin{array}{cc} {{s_{x,PR}}\left( {\frac{{PR}}{2}} \right)}\\ {{s_{y,PR}}\left( {\frac{{PR}}{2}} \right)}\\ {{u_{x,PR}}\left( {\frac{{PR}}{2}} \right)}\\ {{u_{y,PR}}\left( {\frac{{PR}}{2}} \right)} \end{array}} \right]\\= \; \left[ {\begin{array}{cc} {{W_{PR}}}&{{W_{PR}}}\\ { - {V_{PR}}}&{{V_{PR}}} \end{array}} \right]{\left[ {\begin{array}{cc} {{e^{ - \frac{{{\lambda_{PR}}{k_0}{L_{PR}}}}{2}}}}&0\\ 0&{{e^{\frac{{{\lambda_{PR}}{k_0}{L_{PR}}}}{2}}}} \end{array}} \right]^{ - 1}}{\left[ {\begin{array}{cc} {{W_{PR}}}&{{W_{PR}}}\\ { - {V_{PR}}}&{{V_{PR}}} \end{array}} \right]^{ - 1}}\left[ {\begin{array}{cc} {{W_{trn}}}&{{W_{trn}}}\\ { - {V_{trn}}}&{{V_{trn}}} \end{array}} \right]\left[ {\begin{array}{c} {{c_{trn}}}\\ 0 \end{array}} \right]\\= \left[ {\begin{array}{cc} {{W_{PR}}}&{{W_{PR}}}\\ { - {V_{PR}}}&{{V_{PR}}} \end{array}} \right]{\left[ {\begin{array}{cc} {{e^{ - \frac{{{\lambda_{PR}}{k_0}{L_{PR}}}}{2}}}}&0\\ 0&{{e^{\frac{{{\lambda_{PR}}{k_0}{L_{PR}}}}{2}}}} \end{array}} \right]^{ - 1}}{\left[ {\begin{array}{cc} {{W_{PR}}}&{{W_{PR}}}\\ { - {V_{PR}}}&{{V_{PR}}} \end{array}} \right]^{ - 1}}\left[ {\begin{array}{cc} {{W_{trn}}}&{{W_{trn}}}\\ { - {V_{trn}}}&{{V_{trn}}} \end{array}} \right]\left[ {\begin{array}{c} {S_{21}^{({global} )}{c_{inc}}}\\ 0 \end{array}} \right]\\= \left[ {\begin{array}{cc} {{W_{PR}}}&{{W_{PR}}}\\ { - {V_{PR}}}&{{V_{PR}}} \end{array}} \right]{\left[ {\begin{array}{cc} {{e^{ - \frac{{{\lambda_{PR}}{k_0}{L_{PR}}}}{2}}}}&0\\ 0&{{e^{\frac{{{\lambda_{PR}}{k_0}{L_{PR}}}}{2}}}} \end{array}} \right]^{ - 1}}{\left[ {\begin{array}{cc} {{W_{PR}}}&{{W_{PR}}}\\ { - {V_{PR}}}&{{V_{PR}}} \end{array}} \right]^{ - 1}}\left[ {\begin{array}{cc} {{W_{trn}}}&{{W_{trn}}}\\ { - {V_{trn}}}&{{V_{trn}}} \end{array}} \right]\left[ {\begin{array}{c} {S_{21}^{({global} )}W_{ref}^{ - 1}\left[ {\begin{array}{c} {s_x^{inc}}\\ {s_y^{inc}} \end{array}} \right]}\\ 0 \end{array}} \right]\\= [{{A_{PR}}} ]\left[ {\begin{array}{*{20}{c}} {S_{21}^{({global} )}W_{ref}^{ - 1}\left[ {\begin{array}{*{20}{c}} {s_x^{inc}}\\ {s_y^{inc}} \end{array}} \right]}\\ 0 \end{array}} \right].\end{array}$}$$
$${s_{z,PR}}({PR/2} )= \frac{{\left( {\widetilde {{K_x}}{u_{y,PR}}\left( {\frac{{PR}}{2}} \right) - \widetilde {{K_y}}{u_{x,PR}}\left( {\frac{{PR}}{2}} \right)} \right)}}{{j{\varepsilon _{PR}}}},$$
$${I_{\frac{{PR}}{2}}}({x,y} )= {\left|{F{T^2}\left\{ {{s_{x,PR}}\left( {\frac{{PR}}{2}} \right)} \right\}} \right|^2} + {\left|{F{T^2}\left\{ {{s_{y,PR}}\left( {\frac{{PR}}{2}} \right)} \right\}} \right|^2} + {\left|{F{T^2}\left\{ {{s_{z,PR}}\left( {\frac{{PR}}{2}} \right)} \right\}} \right|^2},$$
$$\begin{aligned}{S^{({global} )}}({MASK} )&= {S^{({ref} )}} \otimes {S^{({mask} )}} \otimes {S^{({PMMA} )}} \otimes {S^{({Al} )}} \otimes {S^{({Si{O_2}} )}} \otimes \ldots \otimes {S^{({PR} )}} \otimes {S^{({trn} )}}\\&= {S^{({ref} )}} \otimes {S^{({mask} )}} \otimes {S^{({TRN} )}},\end{aligned}$$
$$\begin{aligned}{\Psi _{PR}}\left( {\frac{{PR}}{2},MASK} \right) &= \left[ {\begin{array}{c} {{s_{x,PR}}\left( {\frac{{PR}}{2},MASK} \right)}\\ {{s_{y,PR}}\left( {\frac{{PR}}{2},MASK} \right)}\\ {{u_{x,PR}}\left( {\frac{{PR}}{2},MASK} \right)}\\ {{u_{y,PR}}\left( {\frac{{PR}}{2},MASK} \right)} \end{array}} \right] = [{{A_{PR}}} ]\left[ {\begin{array}{c} {S_{21}^{({global} )}({MASK} )W_{ref}^{ - 1}\left[ {\begin{array}{c} {s_x^{inc}}\\ {s_y^{inc}} \end{array}} \right]}\\ 0 \end{array}} \right]\\&= [{{A_{PR}}} ]\left[ {\begin{array}{c} {{{[{{S^{({ref} )}} \otimes {S^{({mask} )}} \otimes {S^{({TRN} )}}} ]}_{21}}W_{ref}^{ - 1}\left[ {\begin{array}{c} {s_x^{inc}}\\ {s_y^{inc}} \end{array}} \right]}\\ 0 \end{array}} \right].\end{aligned}$$
$$RMSE = \sqrt {\mathop \sum \limits_{x,y} \frac{{error{{({x,y} )}^2}}}{{{N_x} \times {N_y}}}} ,\textrm{e}rror({x,y} )= \,{I_{M - RCWA}}({x,y} )- {I_{ANSYS}}({x,y} ). $$
$$\begin{array}{l}{k_x} = m\frac{{2\pi }}{{{D_x}}}\;\;\;m = - \infty , \ldots \; , - 2, - 1,\; 0,\; 1,\; 2, \ldots \; \; ,\; \infty,\\{k_y} = n\frac{{2\pi }}{{{D_y}}},\;\;\;n = - \infty , \ldots \; , - 2, - 1,\; 0,\; 1,\; 2, \ldots \; \; ,\; \infty,\end{array}$$
Select as filters
Select Topics Cancel
Top
       
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.
Privacy | Terms of Use