EUV mask model based on modified Born series

: Mask model is a critical part of computational lithography (CL). Owing to the significant 3D mask effects, it is challenging to accurately and efficiently calculate the near field of extreme ultraviolet (EUV) masks with complex patterns. Therefore, a method based on the modified Born series (MBS) was introduced for EUV mask modeling. With comparable accuracy, the MBS method was two orders of magnitude faster than the finite-difference time-domain method for the investigated examples. Furthermore, the time required for MBS was further reduced when the mask pattern was slightly changed. The proposed method shows great potential for constructing an accurate 3D mask model in EUV CL with high efficiency.


Introduction
Extreme ultraviolet lithography (EUVL) has been a promising technology for maintaining Moore's law for decades [1,2].The resolution and quality of EUVL can be significantly improved by computational lithography (CL) techniques, such as optical proximity correction (OPC) and inverse lithography technology (ILT) [3][4][5].OPC utilizes the mask model, vector imaging, resist model, and etch model in a cascade to simulate the lithographic process and predict the achieved wafer patterns, which are modified iteratively until the predicted wafer patterns meet the requirements.Because the mask model is invoked in every iteration, it is one of the most important parts of the CL that determines the performance and efficiency of the OPC.
As the lithography process transcends the advanced nodes, the critical dimension continues to shrink to the scale of the mask thickness, and the 3D mask effect must be considered in the CL.EUVL uses reflective optics, which results from the high light absorbance of almost all materials at an operating wavelength of 13.5 nm.EUV masks were used for reflection under off-axis illumination.Moreover, the height of the absorber was approximately 50 nm, which is greater than the wavelength.These characteristics make EUVL prone to more significant 3D mask effects, such as the asymmetric shadowing effect, deformation of the wavefront, and image contrast variation with feature orientations [6,7].Therefore, it is important to consider 3D mask effects in the EUV mask model.
EUV Mask models have been extensively studied, and can be classified into rigorous and decomposition methods.Rigorous methods, such as the finite-difference time-domain (FDTD) [8], pseudo-spectral time-domain (PSTD) [9,10], finite element method (FEM) [11], boundary element method (BEM) [12], waveguide method (WG) [13], and rigorously coupled wave analysis (RCWA) [14], can provide accurate near-field results for EUV masks.However, their computational cost and time consumption render them unsuitable for large-scale EUV mask simulations.To meet the computational speed requirements, decomposition methods have been proposed to balance accuracy and speed.Domain decomposition method (DDM) decomposes the pattern into simple features and stores the results of these features as a library.Near fields are obtained by synthesizing the precomputed results [15][16][17][18], which limits the accuracy of the complex curvilinear pattern.The structure decomposition method (SDM) [19][20][21][22] models absorbers and multilayers separately.To simulate the absorber, compact models based on the Kirchhoff model with calibrated modifications were proposed.For example, pulses [20] and boundary layers [23] are added to the edge of the absorber.A better accuracy can be achieved using different pulses at different incident angles [22].These compact models are fast but unsuitable for complex patterns.The available methods for the simulation of multilayers are the single-surface approximation [24], ray-tracing method [25] and transfer matrix method (TMM).Simulations of the absorber and multilayer were linked using Fourier transform.
The Born Series is a solution to the Helmholtz equation for the electric field that starts with the famous Born approximation and iteratively calculates the field [26].The Born approximation is limited to weak scattering problems, in which the variation in the refractive index is small.Therefore, the Born series is used in biological applications, such as 3D phase microscopy [27] and optical diffraction tomography [28].Because the Born series is not guaranteed to converge, Osnabrugge et al. proposed a modified Born series (MBS) method by applying a preconditioner [29].Furthermore, Krüger et al. introduced the vector MBS to solve the full wave equation for isotropic situations, which showed an accuracy comparable to that of FDTD [30].
Inspired by the advantages of MBS in solving the weak scattering problem, which is an important characteristic in EUV mask modeling, we introduce the MBS method for EUV mask modeling.Furthermore, because MBS solves Maxwell's equations in an incremental iterative scheme, the efficiency of OPC can be significantly enhanced by employing an MBS-based mask model in which only small changes are applied to the mask patterns in each step.The previous mask pattern can be used as a good initial condition for the near-field calculation of the changed pattern, leading to fewer iterations.Additionally, the proposed method can easily adapt to the requirements of the next-generation EUVL employing a high numerical aperture [1,3] and curvilinear mask [31], which is expected to be a potentially powerful method for EUV mask modeling with a low calculation cost.Detailed descriptions of the model are provided in Section 2. The accuracy and efficiency of the proposed model are discussed in Section 3.

Frame of mask model
A schematic diagram of the EUV mask structure is shown in Fig. 1(a).The EUV masks consisted of a TaN absorber and a Mo/Si multilayer.Because of multilayer reflection, the diffraction of the absorber was calculated twice.The proposed model employs the SDM, and the computational process is illustrated in Fig. 1(b).E in (x,y), E down (x,y), E up (x,y), and E out (x,y) represent the incident, downward, upward, and output fields, respectively.The 3D fields in the downward and upward diffractions are denoted as E diff down (x, y, z) and E diff up (x, y, z), respectively.First, the downward diffraction was simulated using MBS.Details of the MBS are provided in Section 2.3.Note that the field reflected directly by the absorber can be obtained from downward diffraction.However, the reflection from the absorber was significantly smaller than that from the multilayers.This part of the reflection was ignored.E diff down (x, y, z) and E diff up (x, y, z) can be considered the initial conditions for the calculation of the next pattern if the pattern is slightly changed, as indicated by the dashed line in Fig. 1(b).The reflection of the multilayer was calculated using TMM.A detailed description of the TMM and the connection between the TMM and MBS is provided in Section 2.2.Finally, MBS was performed again to simulate the upward diffraction.The described framework of the mask model was used throughout this study.

Transfer matrix method for multilayer
The TMM is a powerful method for the simulation of homogeneous multilayer stacks.The mask model in this study focuses on EUV masks without buried multilayer defects.Theoretically, a multilayer comprises 40 bilayer repetitions.However, the measured reflectance curves of the multilayer did not match the simulated curves because of intermixing between the Si and Mo layers.Therefore, a multilayer configuration with intermixing layers [32] was used in this study.
The TMM was used to calculate the complex reflection coefficients under different incident angles.
The calculation of the multilayer reflection can be divided into five steps.The first step is decomposition. Ẽdown where F represents the forward Fourier transforms; α m and β m represent the direction cosine of the m-th order; the subscripts indicate the components of fields; the tilde represents the Fourier transform of corresponding fields.The second step is the coordinate transformation.The (x, y, z) coordinate system can be transformed into the (s, p) coordinate system [33].The (s, p) system is a local intrinsic system that is used to describe plane waves.
In the third step, the plane-wave components are multiplied by the complex reflection coefficients.
where r s m and r p m represent the s-polarized and p-polarized reflection coefficients of the m-th order, respectively.Subsequently, the (s, p) coordinate system was transformed back into the (x, y, z) coordinate system.
The final step is composition.The processed plane-wave components composes E up (x,y).
where F −1 represents the inverse Fourier transforms.

Modified Born series for absorber
Assuming that the permeability is invariant in space and the field is time harmonic, an integral equation of the scattering potential can be derived from Maxwell's equations by volume integral methods using the Green's function.The MBS is based on an integral equation for potential scattering.A thorough exposure was observed in [26,30,34].The integral equation for the scattering potential is where E(r) is the field at position r ∈ R D ; D is the dimension of the system; g 0 (r) represents the Green's function; S(r) is the source term; V(r) is the scattering potential and V(r is the wavenumber at position r; k 0 is a constant background wave number; ϵ is a constant real number.Equation ( 6) can be spatially discretized and written in matrix form E = GVE + GS (7) where G denotes the convolution in Eq. ( 6) and V is a diagonal matrix containing the V(r).E and S are vectors that contain the values for every cell.The diffraction of the absorber is calculated by solving Eq. ( 7).Practically, it is difficult to solve it in a closed form.In the implementation, the integral equation is solved based on the classical fixed-point iteration [35] where E 0 denotes the initial condition of the iteration.With good initial conditions close to the solution, the number of iterations can be significantly reduced.When E 0 = 0, the iterations are mathematically equivalent to power expansion Equation ( 9) represents the Born series.By truncating the series, we obtained a numerical solution of the integral equation.The Born series converges only if the magnitudes of all eigenvalues of GV are less than one, whereas the MBS converges unconditionally.MBS can be expressed as where γ = (i/ε)V and M = γGV − γ + I.By default, the MBS used in this study refers to a vector-modified Born series.The vector-modified Born series can simulate the polarization effect, which is not negligible for EUV lithography below sub−7 nm generations [36].

Results and discussion
Given that FDTD is a standard method for the simulation of thick mask, a comparison between FDTD and MBS on application of EUV mask absorber is necessary.By default, the discretization in the FDTD and MBS was 28 and 6 points per wavelength (PPW), respectively.In FDTD simulations, a perfectly matched layer was used to absorb the outgoing light.The counterpart of the MBS is the polynomial boundary layer introduced by Osnabrugge In the other directions of the MBS, an antireflection boundary layer [37] with a thickness of two wavelengths at each side was used.The iteration in the MBS stops when the maximum modification of the field between two iterations over all the cells is less than 0.5% of the amplitude of the incident field.These parameters provide a good trade-off between accuracy, simulation time, and memory usage.By default, the initial condition of the MBS is zero.
The absorber material was TaN with n = 0.9385 and k = 0.03776.The height of the absorber layer was 50 nm.The near fields were extracted at the plane, which was 13.5 nm above the absorber.All simulations in this study used a multilayer configuration from Philipsen [32].In most situations, a plane wave is applied to a patterned mask with a nonperiodic boundary.However, applying a plane wave in a nonperiodic domain is possible only if the domain is several orders of magnitude larger than the wavelength, which results from the diffraction of light close to the boundaries caused by the limited simulation domain.A windowed plane wave was used to alleviate the diffraction of light close to the boundaries.Windowed plane waves are generated by applying amplitude modulation to ordinary plane waves [37].The s-polarized light hit the mask at an incident angle of 6°.The transformation between the (s, p) and (x, y, z) coordinate systems is expressed by Eq. ( 4).Only the component E y of the incident fields is not zero, and component E y is the largest component in the near field.Therefore, only the component E y of the near field was demonstrated by default.In addition, both MBS and FDTD represent fields as complex numbers rather than as separated amplitudes and phases, indicating that the amplitudes and phases are treated equally.This was sufficient for comparing the amplitudes of the near fields.
The workstation used for the computation had an Intel Xeon Gold 6285R processor with 56 CPU cores and 512 GB of memory.An Nvidia RTX 5000 with 16 GB of memory is also available.The FDTD solver used in this study is a Lumerical [38] 3D Electromagnetic Simulator.The MBS was implemented using MATLAB.

Accuracy
The calculated pattern is the logical pattern shown in Fig. 2(a).The total simulated area was 810 nm × 810 nm.The dimensions were specified on mask scale by default.The minimum width of the lines in the pattern is approximately 35 nm.The vertical bar on the left is 78 nm wide and has rounded ends to investigate the viability of the MBS on curvilinear EUV masks.Theoretically, the complexity of the curvilinear patterns of EUV masks that can be computed by MBS is limited by discretization.
The convergence of the proposed model is investigated using this pattern.As the PPW increases, the result of the MBS should converge, and the relative error should reach a stable value.The relative error is defined as As shown in Fig. 3, the relative error decreases rapidly at a low PPW and gradually decreases at a high PPW.The convergence of the MBS on discretization was proved.
The proposed mask model using the MBS was compared with the same model using FDTD.The near fields from the MBS and FDTD are shown in Fig. 2  An evident asymmetric effect was observed on the right side of the pattern.The amplitude changes at the edge of the simulation region resulted from windowed plane waves.Figure 2(d) shows that most of the amplitude difference occurred at the edge of the pattern, and the largest amplitude difference is −0.0546.Besides the amplitude, the phase of near fields is also important.The phase difference of near fields calculated by MBS and FDTD is shown in Fig. 4. Phase differences in most of region were small, while phase differences at the positions with small amplitudes were very large because the phases at these locations were sensitive to numerical error.Considering an error of −0.1i, the resulting phase change is 2.7°for a complex number of 0.1 + 0.5i.However, for a complex number like 0.01 + 0.05i, which has a small amplitude, the resulting phase change is 157.4°.Although these phase differences were large, the corresponding errors were relatively small.For clarity of presentation, the range of phase difference was set between −30°and 30°.Actually, both MBS and FDTD represent fields as complex numbers rather than as separated amplitudes and phases, indicating that the amplitudes and phases are treated equally.Only the amplitudes were compared in other examples for brevity.Edge placement error (EPE) resulting from the near field difference was investigated.Usually, EPE is the difference between the intended and the printed pattern, while EPE in this study is defined as the difference between predicted patterns based on near fields from MBS and FDTD.The aerial images, the intensity distribution of field at wafer plane, were calculated with vector imaging [33].The Numerical aperture and demagnification of the simulated aberration-free project system was 0.33 and 4, respectively.Due to the diffraction-limited feature of the projection system, high frequency part of the near field was discarded.As shown in Fig. 5, the aerial images from MBS and FDTD were almost identical.The maximum difference between them was −0.0188, which located at y = 25.88 nm.The intensity of field at y = 25.88 nm is shown in Fig. 5(d).The resist profiles were obtained by a constant threshold resist model with the threshold of 0.5.The EPE at x = 24.75nm and 54 nm was 0.56 nm, which equaled the length of one pixel in the simulation.The EPE at x = −22.5 nm and −60.2 nm was 0 nm.The analyses provide an estimation of EPE resulting from the near field difference.
Another example with a square array pattern is shown in Fig. 6.Compared with the case shown in Fig. 2, only the patterns changed.The width of the square in Fig. 6(a) was 50 nm.The pitch of the square array was 150 nm.The largest difference in this example was 0.0347.Most differences were located at the edge of the pattern.A periodic example with line/space pattern is shown in Fig. 7.The pitch and the width of the absorber were 516.75 nm and 101.25 nm, respectively.The incident field was s-polarized with the incident angle of 6°.The amplitude of the field from MBS and FDTD are shown in Fig. 7, where the dashed line represents the cross-sectional contour of the absorber.The near fields from the MBS and FDTD were consistent at most positions.
As the illumination in EUV lithography is partially coherent, it is necessary to calculate the near field under different source points.The proposed method can be applied to EUV mask with complex curvilinear pattern under various incident angles, without need of calibration or modification.These advantages are demonstrated by example shown in Fig. 8.The total simulated area was 1620 nm × 1282.5 nm.The minimum width of the pattern was 50 nm.The incident field was similar to the one in Fig. 2, but the incident angle was changed to 9°instead of 6°.The maximum amplitude difference was −0.047, which showed great flexibility of the proposed method.Considering the marginal disparities in all instances, MBS emerged as a viable substitute for FDTD in EUV mask absorber simulations.The time consumptions of MBS and FDTD are listed in Table 1.For the logic pattern shown in Fig. 2(a), the MBS simulations are approximately 118 times faster than the FDTD simulations on the CPU.The memory usage of MBS is considerably smaller than that of FDTD, which enables GPU acceleration.An MBS on a GPU is approximately 13 times faster than that on a CPU.In this case, the MBS can calculate three orders of magnitude faster than the FDTD method.Similar results were observed for the other two cases.
The speed advantage of MBS over FDTD in EUV mask simulations results from fewer iterations and coarser discretization.Typically, the iterations of the MBS are small in situations where the variation in the refractive index is small [30].The refractive indices of common materials for EUV mask absorbers are close to those for vacuum, as listed in Table 2, which reduces the number of iterations of the MBS.Additionally, because the discretization of FDTD is constrained by second-order accuracy [8], fine discretization is required for an acceptable error.
In contrast, MBS implements the integrand by a Fourier transform, which has spectral accuracy [39] and allows coarser discretization without sacrificing accuracy.

Expected performance improvement in OPC
The accuracy of the MBS with the initial condition from a mask with a similar pattern was investigated.The modification of the pattern shown in Fig. 9(a) was used to mimic the pattern optimization of the OPC.The pattern was modified by adding hammerheads.A hammerhead is a square pattern added to the outer corner of the original pattern, which is used to compensate for the line-end shortening problem in lithography [41].The width of the hammer head was 9 nm.The amplitudes of the near fields of the modified pattern calculated using FDTD are shown in Fig. 9(b).The difference between the results from FDTD and MBS with the initial condition from result of the original pattern is shown in Fig. 9(c).The difference between the results of FDTD and MBS with zero initial conditions is shown in Fig. 9(d).The differences between the two cases are approximately the same, indicating that the accuracy of the MBS is independent of the initial conditions.Three typical pattern corrections of the OPC, including hammerhead, mousebit, and subresolution assist features (SRAF), were applied individually to demonstrate the performance improvement by the initial condition from the mask with a similar pattern.Mousebits are square patterns subtracted at the inner corners of the original pattern.SRAF are small bars placed close to the primary pattern, which make the isolated pattern features behave lithographically more like a dense pattern.The modified pattern is shown in Fig. 10(a), where the width of the SRAF bar is 10 nm.The iteration reduction for the SRAF with varying lengths is shown in Fig. 10(b).Even when the length of the added SRAF exceeded twice the wavelength (specifically, 30 nm), the reduction in iterations reached 24%.However, the iteration reductions for the hammerhead and mousebit with a width of 9 nm were 40% and 50%, respectively, as shown in Fig. 10(d) and  (d).For larger widths, the iteration reduction converged to a small value of approximately 20%.

Conclusion
In summary, we proposed a fast and accurate mask model for EUV CL.The model employs the SDM, in which the MBS is used to calculate the diffraction of the absorber, and the TMM is used to simulate the reflection of the multilayer.The results of MBS and FDTD are compared on logic, square array, line/space, and complex curvilinear patterns, which prove that MBS is a suitable alternative to FDTD in the application of the EUV absorber.For the investigated examples, a speed improvement of MBS versus FDTD by two orders of magnitude was observed and further speed advantage can be achieved on a GPU.In addition, the iterations of the MBS can be reduced by the initial condition from the field of the mask with a similar pattern.Considerable efficiency improvement of the MBS by the initial condition is demonstrated in three typical OPC situations, including hammerhead, mousebit, and SRAF, with various degrees of change in pattern, which is advantageous to OPC.
Compared with other mask models, the proposed model can simulate EUV masks with complex curvilinear patterns, and exhibits a good balance between speed and accuracy.The proposed mask model makes wide application of CL in EUVL closer.Although the results suggest several advantages, it is crucial to address some potential limitations.For example, MBS based on Fourier transform is difficult to implement parallelly for large-scale simulations.Besides, the discretization is too coarse to accurately represent the thickness.Further research is needed to overcome these limitations.

PatternFig. 1 .
Fig. 1.(a) Schematic diagram of structure of EUV masks.(b) Computation process of the proposed model.The input of MBS consists of the incident field, pattern information and initial condition.The output of MBS is the diffraction field.It takes three steps to calculate the near fields of mask.MBS is used to calculate the downward diffraction.Then the reflection of multilayer is simulated by TMM.Finally, the near fields are obtained from the result of the upward diffraction.
[29] et al.The thickness of the boundary on each side was four wavelengths in the z direction with N = 5 and max |k 2 (r) − k 2 0 | = k 2 0 .
(b) and (d), and are almost identical.

Fig. 2 .
Fig. 2. (a) Logic pattern of EUV mask.Amplitude of the near fields calculated by (b) MBS and (c) FDTD.(d) The difference between (b) and (c).The amplitude of the near fields on (e) x = 0 nm and (f) y = 0 nm.The source is a windowed y-polarized plane wave with the incident angle of 6°.The minimum width of lines in the pattern is about 35 nm.The absorber consists of a 50-nm-thick TaN layer with n = 0.9385 and k = 0.03776.

Figure 2 (
e) and (f) show the amplitudes of the near fields at x = 0 nm and y = 0 nm.The two lines overlap at most positions.

Fig. 4 .
Fig. 4. Phase of the near fields calculated by (a) MBS and (b) FDTD.(c) The difference between (a) and (b).The calculated pattern is the logical pattern shown in Fig. 2(a).

Fig. 5 .
Fig. 5.The intensity of field at wafer plane based on the near fields from (a) MBS and (b) FDTD.(c) The difference between (a) and (b).(d) The intensity and resist profile on y = 25.88 nm.The threshold of the resist model was 0.5.

Fig. 6 .
Fig. 6.(a) Square array pattern of EUV mask.Amplitude of the near fields calculated by (b) MBS and (c) FDTD.(d) The difference between (b) and (c).The width of square is 50 nm.The pitch of the square array is 150 nm for both directions.Other configurations are the same as in Fig. 2.

Fig. 7 .
Fig. 7. Amplitude of near fields for line/space pattern from MBS and FDTD.The black dashed line is the cross-section contour of the absorber.The pitch and the width of absorber were 516.75 nm and 101.25 nm, respectively.

Fig. 8 .
Fig. 8. (a) A complex curlinear pattern of EUV mask.Amplitude of the near fields calculated by (b) MBS and (c) FDTD.(d) The difference between (b) and (c).The incident angle was 9°.

Fig. 9 .
Fig. 9. (a) Change of pattern used to mimic the pattern optimization of OPC.(b) Amplitude of the near fields of the changed pattern calculated by FDTD.(c) The difference in results between FDTD and MBS with initial condition from result of original pattern.(d) The difference in results between FDTD and MBS with zeros initial condition.

Fig. 10 .
Fig. 10.(a) Change of pattern used to mimic the pattern optimization of OPC.The iteration reduction for (b) SRAF, (c) hammerhead and (d) mousebit.The width of SRAF bar is 10 nm.