The ever smaller mask feature sizes with increasingly large designs of semiconductor industry call for the development of next-generation lithography simulators and increasingly aggressive optical resolution enhancement techniques (RETs) ^[1] including modified illumination schemes ^[2] and optical proximity correction (OPC) ^{[3, 4]}. Source mask optimization (SMO) with inverse lithography technology (ILT) ^[5] generates an ILT-optimized mask for each designated illumination with the best combination of source and mask to compensate for the optical diffraction error introduced by the band-limited imaging system, enabling the continuation of the current immersion lithography.

The goal of the SMO is to transfer designed circuit patterns onto silicon wafers, whose performance is heavily dependent on the accuracy of image formation model and the efficiency of the inverse imaging calculation. Over the years, the evolvement of the lithography systems has witnessed the extension from scalar imaging models which are only accurate for systems with low numerical apertures (NAs) less than \(0.4\) ^[6] to vector imaging models. The development of vector imaging can be tracked back from the long celebrated classics by Born and Wolf ^[7], through the revised format of the Hopkins formulation by Yeung ^[6], Flagello ^[8], Pistor ^[9] and the resist vector model by Adam et al. ^[10] to the consistent high-NA optics by Peng ^[11]. Together with the development of illuminations source representations from the arc decomposition ^[12], the meshpoint representation ^[13] to pixelated sources by customized diffractive optical element (DOE) ^[14], various SMO approaches based on scalar imaging model including gradient methods ^[15-17], augmented Lagrangian ^[18] and alike, are given impetus to include vector imaging formation tracking the polarization state as light propagates through the optical components of the projection instruments ^[17]. Meanwhile, level-set approaches have been actively explored in ^{[19, 20]} solving inverse lithography problems, and systematic level-set implementations with appropriate finite-difference schemes are presented in scalar ^[21-23] and vector imaging systems ^{[24, 25]} for photomask synthesis problems. However, the coupling images in these frameworks are usually assumed to be aerial images, while in optical projection lithography intensity distribution within the photoresist or the resist images should be used to account for the transmission and reflection through the dielectric film stack.

In this paper, we present a variational level-set formulation for the source and mask optimization incorporating a regularization energy term preserving the signed distance property combined with an external energy seeking the minimization of the properly designed cost functions to improve pattern fidelity, securing stable level-surface evolution and computation accuracy and eliminating the need of reinitialization. The underlying vector imaging formation is combined with a stratified media model to account for the vector nature of electromagnetic field propagation in the lithography system and to track the reflection from and transmission through the wafer stack. The level-set formulation is reduced to a time-dependent model which is subsequently solved with Polak-Ribiệre-Polyak (PRP) conjugate gradient (CG) method to improve convergence performance. Electrical field caching technique (EFCT) and fast fourier transform (FFT) are employed to significantly reduce the computation load.

With \(\mathrm{N}\) and \({\mathrm{N}}_{\mathrm{s}}\mathrm{ }\)being the sizes of the discretized mask and source matrices, given a target pattern \({\mathbf{I}}_{0}\in {\mathcal{R}}^{\mathrm{N}×\mathrm{N}}\), the goal of the SMO is to find the optimal source \(\hat{\mathbf{J}}\in {{R}^{{{N}_{S}}\times {{N}_{S}}}}\) and mask pattern \(\hat{\mathbf{M}}\in {{R}^{N\times N}}\), which minimizes a designed distance between \(\mathcal{T}\left\{\bullet \right\}\) and \({\mathbf{I}}_{0}\), namely

\(\left( \hat{\mathbf{J}},\hat{\mathbf{M}} \right)=\underset{\hat{J}\in {{R}^{{{N}_{S}}\times {{N}_{S}}}},\hat{M}\in {{R}^{N\times N}}}{\mathop{\min }}\,d\left\{ {{I}_{0}},\mathcal{T}\left\{ \mathbf{J},\mathbf{M} \right\} \right\}\) , (1)

with \(\mathrm{d}\left\{\bullet ,\bullet \right\}\) being an appropriately defined distance metric and \(\mathcal{T}\left\{\bullet \right\}\) being the forward imaging formation.

It is observed that in Equations (9), (16) and (17), the complexity of the convolution operation is superlinear with respect to \(\mathrm{N}\), which is much more computational expensive than the entry-by-entry multiplication. Therefore the major computation load in updating the mask pattern is invoked by the convolution operations. The first strategy to ease the computation load is referred to as the electric-field caching technique (EFCT) ^{[17, 25]}, where \({\mathbf{E}}_{\mathrm{p}}\left({\mathrm{\alpha }}_{\mathrm{s}},{\mathrm{\beta }}_{\mathrm{s}}\right)\) is computed only once in each iteration for the computation of the aerial image \({\mathbf{I}}_{\mathrm{a}}\) and is cached in the memory for the later computation of \({\mathcal{v}}_{\mathbf{i}}\left(\mathbf{r},\mathrm{ }\mathrm{t}\right)\), hence a tradeoff of extra memory to save \({\mathrm{N}}_{\mathrm{s}}×\mathrm{N}×\mathrm{N}\) matrices is made for fast computation. The computation intensity can be further reduced by applying the Fast Fourier Transfrom (FFT) to remove the convolution operations by implementing them in the frequency domain where \({\mathbf{E}}_{\mathrm{p}}\) is computed as \({\mathcal{F}}^{-1}\left[\mathcal{F}\left({\mathbf{H}}_{\mathrm{p}}\right)\mathcal{F}\left(\mathbf{B}⨀\mathbf{M}\right)\right]\) in Equations (9) and (16), while Equation (17) is rewritten as

In this paper, we update mask by employing the widely used Polak-Ribiệre-Polyak (PRP) conjugate gradient (CG) ^[34] which is believed most efficient, instead of the steepest descent (SD) practiced in Reference [21, 23, 24] to achieve faster convergence and improved performance, where the evolution velocity \(\mathrm{G}\left(\mathbf{r},\mathrm{ }{\mathrm{t}}_{\mathrm{k}}\right)\) in the \(\mathrm{k}\)th iteration is calculated according to Equation (18) as

\(G\left( r,{{t}_{k}} \right)=\left\{ \begin{matrix} g\left( r,{{t}_{k}} \right)+\eta _{k}^{PRP}\cdot G\left( r,{{t}_{k-1}} \right) \\ g\left( r,{{t}_{k}} \right) \\\end{matrix} \right.\begin{matrix} if\ \ k\ge 1 \\ if\ \ k=0 \\\end{matrix},\) (20)

with \({\mathrm{\eta }}_{\mathrm{k}}^{\mathrm{P}\mathrm{R}\mathrm{P}}\) calculated as

\({\mathrm{\eta }}_{\mathrm{k}}^{\mathrm{P}\mathrm{R}\mathrm{P}}=\frac{{\mathrm{g}\left(\mathbf{r},\mathrm{ }{\mathrm{t}}_{\mathrm{k}}\right)}^{2}-\sum \mathrm{g}\left(\mathbf{r},\mathrm{ }{\mathrm{t}}_{\mathrm{k}}\right)\bullet \mathrm{g}\left(\mathbf{r},\mathrm{ }{\mathrm{t}}_{\mathrm{k}-1}\right)}{{\mathrm{g}\left(\mathbf{r},\mathrm{ }{\mathrm{t}}_{\mathrm{k}-1}\right)}^{2}},\) (21)

The update of the time-dependent \({\Phi }\left(\mathbf{r},\mathrm{t}\right)\) can be implemented by available first-order temporal accurate and second-order accurate spatial finite-difference schemes ^[21]. With \({\Phi }\left(\mathbf{r},\mathrm{t}\right)\) discretized as \({{\Phi }}_{\mathrm{i},\mathrm{j}}^{\mathrm{k}}\) using spatial index \(\left(\mathrm{i},\mathrm{j}\right)=1,\mathrm{ }\cdots ,\mathrm{ }\mathrm{N}\) and temporal index \(\mathrm{k}\), the level-set equation in Equation (13) can be expressed as

\({{\Phi }}_{\mathrm{i},\mathrm{j}}^{\mathrm{k}+1}={{\Phi }}_{\mathrm{i},\mathrm{j}}^{\mathrm{k}}+\mathrm{\delta }\mathrm{t}\bullet L\left({{\Phi }}_{\mathrm{i},\mathrm{j}}^{\mathrm{k}}\right), k=0, 1, 2, \cdots , \) (22)

which is a iteration process with \(L\left({{\Phi }}_{\mathrm{i},\mathrm{j}}^{\mathrm{k}}\right)\) being a discretized approximation of \({-\left|\nabla {{\Phi }}_{\mathbf{i}}\right|\mathrm{g}}_{\mathbf{i}}\left(\mathbf{r},\mathrm{ }\mathrm{t}\right)\) in Equation (18) and\(\mathrm{ }\mathrm{\delta }\mathrm{t}\) being the time step. The updating procedure of level-set based PRP-CG method for the proposed SMO approach is described in Table 1

.

Simulations are performed on a partially coherent imaging system with an annular illumination source \({\mathbf{J}}_{0}\) whose outer radius is \({\mathrm{\sigma }}_{\mathrm{o}\mathrm{u}\mathrm{t}}=0.9\) and inner radius is \({\mathrm{\sigma }}_{\mathrm{i}\mathrm{n}}=0.6\), given in Figure 3(a), with the desired target pattern \({\mathbf{I}}_{0}\) in Figure 3(b). The illuminating wavelength is set at \(193\mathrm{n}\mathrm{m}\) and the numerical aperture is \(1.35\). The resolution of the mask pattern is \(\mathrm{\delta }\mathrm{x}=\mathrm{\delta }\mathrm{y}=4\mathrm{n}\mathrm{m}\). The resist effect is approximated by a sigmoid function with steepness \(\mathrm{a}=85\) and the threshold \({\mathrm{t}}_{\mathrm{t}}=0.3\). The parameters of the wafer stack is referenced from ^[2] and is given in Table 2. The multiplier \(\mathrm{\mu }=0.1\) is used in the simulations. All the computations are carried out using MATLAB on an Intel(R) Core(TM) i\(7\)-\(6700\) CPU, \(3.40\) GHz, \(8.00\) GB of RAM and stop after \(150\) iterations. In Figure 3, lithographic image formations with coupling image as the aerial image and the photoresist image are depicted. Figure 3(c) and (e) show the coupling image using vector imaging without and with the stratified model to include the reflection from and transmission through the wafer stack \(10\mathrm{n}\mathrm{m}\) below the interface between the immersion liquid and the wafer stack surface where \(\mathrm{z}=0\) in Figure 2, respectively. Figure 3(d) and (f) give the corresponding wafer images of Figure 3(c) and (e). It is observed that severe distortions invoked by the low pass nature of the pupil function \(\mathbf{h}\), exhibiting pattern errors (PEs) \(2361\) and \(8575\) in the simulations. In Figure 4, the same mask pattern in Figure 3(b) is illuminated by a dipole source in Figure 4(a) to form the aerial image in Figure 4(b) with the corresponding wafer image of PE \(2742\) in Figure 4(b), and the coupling photoresist image in Figure 4(d) with the corresponding wafer image of PE \(8343\) in Figure 4(e). Pattern fidelities in Figure 3(f) and Figure 4(e) are further deteriorated because of the intensity absorption of electromagnetic filed propagation in the wafer stack, which has to be compensated by more radical computational techniques.

Figure 5 presents the optical proximity corrections by applying the proposed SMO approach which starts the updating process with different initial source patterns. Row (a) of Figure 5 shows the lithographic imaging performance with the starting annular source \({\mathbf{J}}_{0}\) to generate the synthesized source pattern \(\hat{{\mathbf{J}}_{0}}\) and mask pattern \(\hat{\mathbf{M}}\) and improves the pattern fidelity from PE \(8575\) to \(910\), while the simulation in row (b) starts the optimization process with the dipole source \({\mathbf{J}}_{1}\) and generates synthesized source pattern \(\hat{\mathbf{J}}_{1}\) and mask pattern \(\hat{\mathbf{M}}\) with improved pattern fidelity from PE \(8343\) to \(409\) in the wafer image. The proximity-corrected wafer images in row (a) and (b) are comparable except for the missing contact areas in row (a). It should be noted that the inverse lithography problem is non-convex with multiple local minima. With the time-dependent model in Equation (18) solved with finite-difference schemes, there is no guarantee of reaching global minimum. However, ILT is an ill-posed problem and a global minimum is often not required. Any good local minimum (where goodness is defined by data fidelity and user-defined properties) will suffice as an acceptable solution.

The proposed SMO approach is further applied to another desired mask pattern \({\mathbf{I}}_{1}\) in Figure 6(a) illuminated by the annular source \({\mathbf{J}}_{0}\) in Figure 3(a), where pattern fidelity degradation can be observed from the wafer image of the aerial image with PE \(4363\) in Figure 6(c) and the wafer image of the coupling image with PE \(6590\) in Figure 6(e). OPC with mask optimization (MO) and the proposed SMO approach are performed in the simulations of row (a) and (b) of Figure 7, respectively, where in MO, only the mask pattern is modified to improve the pattern fidelity from PE \(6590\) to \(1280\), while in SMO, both the source and mask patterns are synthesized to further reduce the PE to \(1068\). The pattern fidelity improvement using the proposed SMO approach is justified by the fact that source and mask optimization (SMO) expands the solution space of the source and mask patterns by the joint optimization of the illumination and mask shapes, whereas mask optimization (MO) procedure has limited degree of freedom with fixed source patterns.

Evaluation of runtime with respect to iteration number for simulations in Figure 7 and conventional steepest descent (SD) method for SMO without the efficiency enhancement techniques in Section 3 is given in Figure 8, where the SD method is terminated in \(80\) iterations because of the huge computation load and linear dependence observed. Pattern error convergence for simulations in Figure 7 and conventional steepest descent (SD) method with respect to runtime is evaluated in Figure 9, exhibiting a much more accelerated convergence for the proposed method than the conventional method. A more quantitative comparison of the runtimes is given in Table 3. It is observed that the average runtime in every iteration using the SD method amounts to \(1.105\) hours, accumulating the total runtime to \(88.37\) hours for the whole \(80\) iterations, which is almost unbearably time-taxing. The average iteration runtimes applying the efficiency enhancement techniques in Section 3 to MO and SMO are \(0.0169\) and \(0.0539\) hours, respectively, adding up to \(2.53\) and \(8.09\) hours for the whole simulations, justifying the fact that SMO approach is updating both the source and mask patterns iteratively instead of only updating the mask pattern in the MO approach. Consequently, a convergence acceleration of about \(20\) times is achieved with the proposed method.

It is worth mentioning that free-form source patterns leading to a gray-scale representation in the synthesized sources \(\stackrel{^}{\mathbf{J}}\) can be readily achieved by applying a parametric transformation to reduce the binary constrained problem to an unconstrained one ^[17] with comparable performances in terms of pattern fidelity improvement to the binary source patterns optimized in the proposed method. Also, the CFL condition of the finite difference schemes states that the Euler time step is confined by the boundedness of \({\mathrm{g}}_{\mathrm{i}}\) in the \(\mathrm{x}\) and \(\mathrm{y}\) direction in Equation (19). Euler time steps \({\mathrm{\delta }\mathrm{t}}_{{\mathrm{\varphi }}_{\mathbf{J}}}\) for updating \(\stackrel{^}{\mathbf{J}}\) and\(\mathrm{ }{\mathrm{\delta }\mathrm{t}}_{{\mathrm{\varphi }}_{\mathbf{M}}}\) for updating \(\stackrel{^}{\mathbf{M}}\) are presented in Figure 10, with a respective average Euler time step of \(0.1346\) and \(0.0054\) for the proposed level-set based SMO approach. It is observed that the Euler time step \({\mathrm{\delta }\mathrm{t}}_{{\mathrm{\varphi }}_{\mathbf{M}}}\) for updating the mask pattern is particularly small which we believe inhibits the emerging assistant features when performing the binarization of the updated pixels in the proposed method.

We have presented a variational level-set formulation for lithographic source and mask optimization incorporating an external energy term minimizing the sum of the mismatches between the printed image and the desired one over all locations and a distance metric preserving the closeness of the evolution function \( {{\Phi }}\) to a signed distance function, which eliminates the necessity of re-initialization in every iteration securing stable curve evolution and ensure desirable results. A more accurate vector imaging formation is established to incorporate the stratified media depicting the reflection from and transmission through of the electromagnetic field propagation in the wafer stack. Electric-field caching technique (EFCT) and fast Fourier transform (FFT) are employed to significantly ease the computation load incurred by the convolutions operations in the optimization process, and the Polak-Ribiệre-Polyak (PRP) conjugate gradient (CG) method to improve convergence. Numerical results demonstrate stable evolution of the level set functions with respect to the source and mask patterns, with improved pattern fidelity and faster convergence. The incorporation of the signed distance regularization term ensures stable LSF evolution and accurate computation without the need of reinitialization. Besides, with each entry of \( {{\Phi }}\) denote the signed distance to the zero contour of the LSF, localizing the level surface evolution in the vicinity of the zero-level set to enhance computation efficiency is potentially available with a simple and straight narrow-band implementation of the variational model, leading our future work to include narrow-band methods with large time steps improving metrics such as EPE, PV band. The theoretical and numerical analysis of the proposed ILT brings algorithmic insights into the exploitation of vector imaging model for SMO in next generation lithography with \(22\mathrm{n}\mathrm{m}\) technology node and beyond.