Accelerating extreme ultraviolet lithography simulation with weakly guiding approximation and source position dependent transmission cross coefficient formula

Abstract. Background Mask three-dimensional (3D) effects distort diffraction amplitudes from extreme ultraviolet masks. In a previous work, we developed a convolutional neural network (CNN) that predicted distorted diffraction amplitudes very fast from input mask patterns. Aim In this work, we reduce both the time for preparing the training data and the time for image intensity integration. Approach We reduce the time for preparing the training data by applying weakly guiding approximation to 3D waveguide model. The model solves Helmholtz type coupled vector wave equations of two polarizations. The approximation decomposes the coupled vector wave equations into two scalar wave equations, reducing the computation time to solve the equations. Regarding the image intensity integration, Abbe’s theory has been used in electromagnetic (EM) simulations. The transmission cross coefficient (TCC) formula is known to be faster than Abbe’s theory, but the TCC formula cannot be applied to source position dependent diffraction amplitudes in EM simulations. We derive source position dependent TCC (STCC) formula starting from Abbe’s theory to reduce the image intensity integration time. Results Weakly guiding approximation reduces the time of EM simulation by a factor of 5, from 50 to 10 min. STCC formula reduces the time of the image intensity integration by a factor of 140, from 10 to 0.07 s. Conclusions The total time of the image intensity prediction for 512  nm×512  nm area on a wafer is ∼0.1  s. A remaining issue is the accuracy of the CNN.


Introduction
High-aspect absorbers used in extreme ultraviolet (EUV) masks induce several mask threedimensional (3D) (M3D) effects such as critical dimension (CD) error and edge placement error. 1,2It is necessary to include M3D effects in EUV lithography simulations.M3D effects are caused by the distorted diffraction amplitude from an EUV mask.4][5] However, these calculations are highly time-consuming, especially for optical proximity correction (OPC) applications.
To speed up the EM simulations, several approximation models such as "domain decomposition method" 6,7 and "M3D filter" 8,9 were proposed, which decomposed a mask pattern into two-dimensional (2D), one-dimensional (1D), and zero-dimensional (0D) patterns.In these models, the EM field of a mask pattern was calculated by superposing the EM fields of 2D, 1D, and 0D patterns.][11] An implicit assumption of these models is that the mask pattern is large and isolated.However, EM interaction is nonlocal.The amplitude is affected by the surrounding patterns.In some approximation models, the first-order cross talks between neighboring edges are included, but higher-order cross talks required in the "rigorous domain decomposition method" 12 are neglected.In the case of OPC masks, the pattern densities are high because the main pattern is decorated by many serifs and assist features.Also, advanced OPC mask patterns are curvilinear.It could be difficult to apply these approximation models to OPC masks.
Recently, many attempts have been made to simulate the M3D effects using deep neural network (DNN) such as convolutional neural network (CNN) or generative adversarial network (GAN).They are classified into three types depending on the target of DNN: the near-field amplitude at the object plane, the image intensity at the image plane, and the far-field amplitude at the pupil plane.
4][15][16] However, the near-field amplitude has local oscillation, which makes it difficult to define the loss function of CNN.Also, since the near-field amplitude depends on the incident angle, these models need many CNNs for different source positions.
The image intensity is a natural target of DNN, and GANs have been applied to reproduce the image intensity from the input mask pattern. 17,18However, in these models, the source shape and the aberrations including the defocus are fixed.In OPC applications, the model needs to be reconstructed when the source shape is changed.
In our previous works, [19][20][21] a CNN model has been developed that predicts the far-field diffraction amplitude from the input mask pattern.Our CNN model can be applied to arbitrary mask patterns.Although CNN prediction time is very short, preparing training data by EM simulation takes a long time.In this work, we apply weakly guiding approximation to 3D waveguide model, 5 one of the EM simulation models, which solves Helmholtz-type coupled vector wave equations.By using the weakly guiding approximation, the coupled vector wave equations are decomposed into two scalar wave equations, reducing the computation time to solve the equations.
The diffraction amplitudes calculated by EM simulations depend on the source position.Hopkins' transmission cross coefficient (TCC) formula, which is conventionally used in optical lithography simulations, cannot handle the source position-dependent diffraction amplitudes. 22herefore, in EUV lithography simulations, Abbe's theory has been used to calculate the image intensity.However, the computation time using the TCC formula is much shorter than the time using Abbe's theory, because TCC can be precalculated before the image intensity integration.Since the diffraction amplitude in our model is described in frequency space, the model can be easily incorporated with the TCC formula.In this work, we derive a source position-dependent TCC (STCC) formula starting from Abbe's theory to reduce the image intensity integration time.
In Sec. 2, we explain the architecture of our CNN.In Sec. 3, we apply the weakly guiding approximation to 3D waveguide model.In Sec. 4, we derive STCC formula.Section 5 is the summary.

CNN for Fast EUV Simulation
0][21] Figure 1 is the schematic view of the diffraction amplitudes Aðl; m; l s ; m s Þ from an EUV mask.We show here the vector potential A. Inside the vacuum the vector potential is converted to the electric field E by the following equation: 19 where k and k represent the wave vector and the wave number, respectively.The diffraction amplitude Aðl; m; l s ; m s Þ is divided into the thin mask amplitude [Fourier transform (FT) of the mask pattern] A FT ðl; mÞ and the M3D amplitude A 3D ðl; m; l s ; m s Þ E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 7 ; 4 6 7 The M3D amplitude for each diffraction order ðl; mÞ smoothly depends on the source position ðl s ; m s Þ as shown in Fig. 2. We assume periodic boundary condition with the mask size L. When the mask pattern is periodic, the momentum (or spatial frequency) of the far-field diffraction amplitude ðk x ; k y Þ has discrete number ðk x ; k y Þ ¼ 2π∕Lðl; mÞ.For the convenience of numerical calculation, we also discretize the source position ðs x ; s y Þ ¼ 2π∕Lðl s ; m s Þ.
We assume the maximum source size σ ¼ 1.The source area σ > 1 corresponds to the darkfield illumination, but we do not use this area in lithography.The source position and the diffraction order are restricted by the source shape and the pupil shape as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 7 ; 2 9 8 where NA = 0.33 is the numerical aperture of the projection optics and λ ¼ 13.5 nm is the wavelength.The magnification of the projection optics is 1/4.Only the overlapping area between the pupil and the source has possibility to contribute to the image intensity.At the center of the overlapping area, ðl s ; m s Þ = ð−l∕2; −m∕2Þ as shown in Fig. 2.
We approximate the M3D amplitude by a linear function of source position as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 4 ; 6 8 8 where a 0 is the average of the amplitude in the overlapping area, and a x and a y are the slopes of the amplitude in x and y directions, respectively.We call these three numbers as M3D parameters.In Eq. ( 5), A 3D x ðl; m; l s ; m s Þ is expanded at the center of the overlapping area ð−l∕2; −m∕2Þ.Therefore, inside the overlapping area, a x ðl; mÞ ðl s þ l∕2Þ+a y ðl; mÞ ðm s þ m∕2Þ is small.This improves the accuracy of STCC formula in Sec. 4.
There is another reason using Eq. ( 5).Mask 3D parameters are derived by least square fitting to the amplitudes at the grid points inside the overlapping area in Fig. 2. The larger ðl; mÞ is, the smaller the number of the grid points inside the overlapping area.If the number of the grid points is too small, the overlapping area becomes a line or just a point.In such case, we approximate the amplitude in the area by using only a 0 ðl; mÞ as the average of the amplitude and do not use a x ðl; mÞ and a y ðl; mÞ.Therefore, a 0 ðl; mÞ should represent the average of the amplitude in the overlapping area.
M3D parameters are determined by the mask pattern.Recently, CNN is widely used as pattern recognition technique.In the previous works, [19][20][21] we constructed a CNN that predicted M3D parameters from an input mask pattern (Fig. 3).
3 Weakly Guiding Approximation of 3D Waveguide Model

Mask Clip Size for High Coherent Illumination
In the previous works, 20,21 we used a periodic mask pattern with 256 nm × 256 nm area on the wafer.We assumed that the area was clipped from large mask data.We should not use the edges of the clipped area to avoid the influence of the neighboring mask pattern.According to Ref. 23, the optical interaction range R opt is calculated by the following equation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 4 ; 3 8 5 where λ; σ, and NA represent the wavelength, coherence factor, and numerical aperture of the scanner, respectively.The wavelength of EUV light is 13.5 nm and the numerical aperture of the current EUV scanner is 0.33.The coherence factor depends on the illumination setting.Equation ( 6) can be used for conventional illumination.Here, we confirm the validity of the equation for high coherent illumination, such as dipole illumination.Figure 4 shows the pitch dependence of 20 nm line CD for conventional and dipole illumination.We use a simple threshold model fixing the threshold intensity value at 40 nm pitch.CD varies depending on the pattern pitch, but it becomes stable at larger pitches where the line is isolated from the neighbor lines.The minimum pitch where CD becomes stable depends on the illumination.The minimum pitch physically corresponds to the optical interaction range.From Fig. 4, the optical interaction range for σ 0.5 is ∼100 nm and that for dipole illumination (outer σ 0.7, inner σ 0.3, and open angle 90 deg) is ∼150 nm.The value for σ 0.5 is close to R opt ∼ 90 nm in Eq. ( 6).In the case of the dipole illumination, the size of each monopole is∼0.3.From Eq. ( 6), R opt for σ 0.3 is ∼150 nm, and it is same as the value derived from Fig. 4.
Fig. 3 CNN that connects a mask pattern and M3D parameters.
Figure 5 shows the usable mask area excluding the area influenced by the neighboring mask pattern.The mask clip size L should be larger than 2 × R opt to get usable mask area.Therefore, when we use high coherent illumination, the mask clip size L needs to be larger than 300 nm.In the previous works, we used L ¼ 256 nm, but it was not large enough for high coherent illumination.In this work, we enlarge the mask clip size to 512 nm to obtain usable mask area for high coherent illumination.The usable area on the wafer is ∼200 nm × 200 nm.

Weakly Guiding Approximation
When we enlarge the mask clip size, the computation time for EM simulations increases.We use 3D waveguide model 5 to solve Maxwell's equations.The calculation time for a 256 nm × 256 nm mask clip is 146 s and the time for a 512 nm × 512 nm mask clip is 2,850 s by using Core i9-10940 central processing unit.The model slices an EUV mask into multilayers including absorber layers and Mo/Si reflective layers.Inside each layer Maxwell's equations are reduced to the Helmholtz type coupled vector wave equations as follows: E Q -A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 7 ; 1 5 3 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 7 ; 1 0 4 where A x and A y are the x and y components of the vector potential.ε is the complex dielectric constant of the absorber layers or reflective layers.Inside each layer, the complex dielectric  constant ε is uniform in the z direction.In this case, gauge transformation freedom allows fixing A z to be zero. 24Then A z is omitted from the coupled vector wave equations, Eqs. ( 7) and ( 8).Inside absorber layers, ε is a function of x and y because the absorber is patterned.Inside reflective layers, ε is uniform in the x and y directions, so Eqs.( 7) and ( 8) can be solved analytically.Two variables, A x and A y , correspond to two polarizations.Equation (1) indicates that the electric fields E of A x and A y polarizations are almost parallel to x and y axes because k x ; k y ≪ k near the optical axis.Figure 6 is an example of diffraction amplitudes calculated by solving Eqs. ( 7) and ( 8) (for the details, see Ref. 19).The result shows that the polarization change between the incident wave and the outgoing wave is very small.This is because the complex dielectric constant of EUV absorber is close to one.Similar phenomenon is known as "weakly guiding approximation" in optical fiber, 25 where two polarizations are decoupled.
We apply the weakly guiding approximation to 3D waveguide model and decompose the coupled vector wave equations.Each equation becomes a scalar wave equation as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 4 ; 5 7 9 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 4 ; 5 3 3 Each equation can be solved independently, and it takes 289 s for a 512 nm × 512 nm mask clip.Solving two equations take 578 s and it is ∼1∕5 of the time for solving original 3D waveguide model.
We confirm the accuracy of the weakly guiding approximation.Figure 7 compares the image intensities calculated by using the 3D waveguide model and the weakly guiding approximation.We assume the conventional Ta absorber, where the complex refractive index ðn; kÞ ¼ ffiffi ffi ε p ¼ ð0.9567; 0.0343Þ. 26It is close to the complex refractive index of the vacuum (1, 0).The difference between the 3D waveguide model and the weakly guiding approximation is very small, <0.1%.Polarization changes due to the EUV mask are negligible.However, there is small difference between A x and A y polarizations.It is expected that the difference becomes larger for high NA scanners where the incident angle becomes large.The result here shows the polarization effect on the mask.We do not include the polarization effect at the exit pupil of the projection optics, which is significant in high NA optics.
Figure 8 shows the results when a low-n absorber TP1 in Ref. 27 is used.The complex refractive index of TP1 absorber is (0.91, 0.032).As shown in Fig. 8, the difference between the 3D waveguide model and the weakly guiding approximation becomes large, at most 1.4%.The accuracy of the weakly guiding approximation is deteriorated when low-n absorbers are used.Note that low-n absorbers are still under development and the mask process has not yet been established as discussed in Ref. 27.
Fig. 6 Polarization dependence of the diffraction amplitudes calculated by 3D waveguide model.

STCC Formula 4.1 Thin Mask Model and Thick Mask Model
According to the Abbe's theory, the image intensity of the thin mask model I Thin is calculated by the following equation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 7 ; 1 9 6 where S is the effective source and P is the pupil function of the projection optics.P is a matrix for high NA optics, 28 but we assume here as a scalar function.The electric field of the thin mask E FT is calculated from the vector potential of the thin mask A FT by using Eq. ( 1).Hopkins' TCC formula 22 is derived by interchanging the order of the integration in Eq. ( 11).
E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 7 ; 1 1 2 where  TCCðk; TCC does not depend on the mask pattern.The calculation time of the image intensity is reduced by precomputing and tabulating TCC.
Abbe's theory is valid in the thick mask model, and the image intensity I Thick is calculated as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 1 1 4 ; 5 3 1 The electric field of the thick mask E is calculated from the vector potential of the thick mask A in Eq. ( 2).The electric field depends on the source position s.Therefore, we cannot interchange the order of the integrations in Eq. ( 14).We cannot apply Hopkins' TCC formula to the thick mask model.Figure 9 compares the image intensities of the thick mask model and the thin mask model.We use the same mask pattern, the same illumination, and the same absorber as shown in Fig. 7.
The maximum difference between the thick mask model and the thin mask model is large, 5.2%.
The shadowing effect at the edges of the absorber is clearly seen in Fig. 9.

Linear Approximation of the Thick Mask Model
According to Eqs. ( 2) and ( 5), the vector potential of the thick mask model is approximated by a linear function of the source position.The electric field of the thick mask model is also approximated by a linear function of the source position as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 1 1 4 ; 3 5 1 Inserting Eq. ( 15) into Eq.( 14), we obtain the image intensity I Linear for the linear approximation of the thick mask model as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 1 1 4 ; 3 0 0 Figure 10 compares the image intensities of the thick mask model and its linear approximation.The maximum difference is 1.4% which is about 1/4 of the difference between the thick mask model and the thin mask model.Linear approximation is a good approximation as a starting point to include the M3D effects, but we might need to include higher order terms if higher accuracy is required.

STCC Formula
The STCC formula that includes the source position dependence of the diffraction amplitudes is derived by interchanging the order of integrals in Eq. ( 16) as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 1 1 7 ; 7 0 0 where TCC x and TCC y are defined by the following equations: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 8 ; 1 1 7 ; 5 6 6 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 1 1 7 ; 5 1 8 In Eq. ( 17), we ignore the second-order term of ∂ s because the contribution is small inside the overlapping area of the source and the pupil, as discussed in Sec. 2.
Sum of coherent systems (SOCS) model 29 is conventionally used in optical lithography simulations to speed up the image intensity integration.SOCS model decomposes TCC into eigen functions and sums up only small number of the eigen modes to calculate the image intensity.SOCS model can also be applied to TCC x and TCC y because they are Hermitian matrices.Then, three TCCs are written as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 1 1 7 ; 4 0 1 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 1 ; 1 1 7 ; 3 5 5 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 2 ; 1 1 7 ; 3 2 6 TCC y ðk; where α n ; β n , and γ n are eigen values, and φ n ; ϕ n , and ψ n are eigen functions.These eigen values are real numbers.Figure 11 compares the image intensities calculated by the linear approximation of the thick mask model using Abbe's theory and STCC formula.In the case of STCC formula, we use SOCS model with 100 eigen modes for TCC, and 20 eigen modes for TCC x and TCC y .The difference between the image intensities calculated by the two formulas is very small, <0.3%.STCC formula reduces the computation time.The computation time of image intensity integration by Abbe's theory is 10 s for 512 × 512 points.On the other hand, the computation by STCC formula takes only 0.07 s excluding the time for the eigen value decomposition by SOCS model.In this work, we accelerated the EUV lithography simulation based on the CNN model.A remaining big issue is the accuracy of the CNN. 21The accuracy depends on the quality and quantity of the training data.We hope that large-scale training mask data will improve the accuracy of the CNN.This work is based on the prior SPIE proceedings paper. 30

Fig. 1
Fig. 1 Schematic view of light diffraction by an EUV mask.Diffraction amplitudes depend on both the diffraction order and the source position.

Fig. 2
Fig. 2 Source position dependence of M3D amplitude.The diffraction order of the amplitude is (l; m).The center of the overlapping area between the source and the pupil is (l s ; m s ) = (−l∕2; −m∕2).

Fig. 8
Fig.8Image intensities calculated by the 3D waveguide model and the weakly guiding approximation.We use a 45 nm thick TP1 absorber.

Fig. 7
Fig. 7 Image intensities calculated by the 3D waveguide model and the weakly guiding approximation.The dipole illumination has σ in ∕σ out ¼ 0.55∕0.9and open angle = 90 deg.We use a 60 nm thick Ta absorber.

Fig. 9 Fig. 10
Fig. 9 Image intensities calculated by the thick mask model and the thin mask model.

Fig. 11
Fig. 11 Image intensities calculated by the linear approximation of the thick mask model (Abbe's theory) and STCC formula.

Figure 12
Figure12shows the estimation of the runtime for CNN preparation, training, and prediction.Weakly guiding approximation reduces the time of EM simulation by a factor of 5, from 50 to 10 min.STCC formula reduces the time of the image intensity integration by a factor of 140, from 10 to 0.07 s.The total time of the image intensity prediction for 512 nm × 512 nm area on wafer (usable area: 200 nm × 200 nm) is ∼0.1 s.In this work, we accelerated the EUV lithography simulation based on the CNN model.A remaining big issue is the accuracy of the CNN.21The accuracy depends on the quality and quantity of the training data.We hope that large-scale training mask data will improve the accuracy of the CNN.This work is based on the prior SPIE proceedings paper.30