Rigorous Electromagnetic Field Modeling

The dimensions of mask features on state-of-the-art lithography masks are small compared to the illumination wavelength. In the EUV case currently mask features are indeed larger but also significantly thicker than the illumination wavelength. Additionally things like sub-wavelength assist features in aggressive OPC, edges on alternating PSM, complicated multilayer structures on EUV-masks, and other geometrical details must be taken into account. Therefore the standard Kirchhoff approach for the mask description cannot be used anymore. A rigorous computation of the light diffracted by the mask must be performed.

In Dr. LiTHO two different approaches for rigorous mask simulations are available: The FDTD (finite-difference time-domain) method and the waveguide method. In general, the diffraction of light from a mask is considered as an electromagnetic boundary value problem. The electric and magnetic fields in the vicinity of the mask are described by Maxwell’s equations. Most materials used in semiconductor lithography are isotropic, non-magnetic and source-free. The basic difference between the methods is the fact that FDTD is operating in the time-space domain and waveguide in the frequency domain.

The FDTD method is widely-used and well established at Fraunhofer IISB since several years. The relatively new implemented waveguide method is a very powerful alternative. At Fraunhofer IISB currently waveguide is used by default for litho simulations. The well established FDTD method serves as reference and is used for special applications.

FDTD method

The basic idea of FDTD is to numerically integrate the first two Maxwell equations in discrete time steps. Both, electric and magnetic field are composed of 3 field components. In other words, 6 mutually dependent field components have to be integrated over time. This integration is performed on equidistant numerical grids. Every field component is specified on an individual grid. Appropriate positioning of the meshes of the different field components guarantees, that the last two Maxwell equations are also fulfilled during the numerical integration procedure.

The accuracy of a numerical solution obtained with FDTD is governed by the spatial discretization dx. As a rule of thumb a relative accuracy better than 2% requires N > 15 meshpoints per wavelength in the material with the largest refractive index n. The stability conditions for the finite difference scheme limits the magnitude of discrete time steps to dt < dx/c, where c is the vacuum velocity of light.

Meshpoints on the boundary require specific updating equations which determine the behavior of the electromagnetic field at the boundaries of the computing window. Possible choices are: Reflecting boundary conditions, periodic boundary conditions, and transparent boundary conditions.

The standard FDTD scheme is only stable for positive permittivities. In other words, the standard FDTD formulation will fail for materials with a strong absorption. The introduction of a special updating scheme for absorbing materials leads to a method which is stable for virtually all materials. But it requires the updating of up to 3 additional field quantities. Therefore, strongly absorbing materials in the computing window require more memory and computing time. For the extraction of steady state fields the control of convergence is either done by comparison with the field values during previous time steps or by certain experience based rules.

Waveguide method

Figure 1 shows the basic idea of the computation of light diffraction from a lithography mask with the waveguide method. For a better understanding the figure shows only a two dimensional cross section of the real three dimensional problem.

Computation of Light Diffraction
Figure 1: Computation of light diffraction from a litho mask using the waveguide method

Based on the real mask geometry a slicing of the area to be simulated is performed by defining maximum sized layers which are homogeneous in z-direction (this is the waveguide assumption, see layers in Figure 1). The geometry (material distribution) of the individual layers is described by an arbitrary number of rectangles with arbitrary sizes (depending on the geometry and on the required resolution). With this mask description a potential sampling error can be limited to any value. Then the material distributions and the electromagnetic fields of all layers are described by Fourier series. All Fourier series are truncated according to the number M of modes which are supposed to be taken into account (number of reflected and transmitted plane waves in Figure 1). The combination of this material and field descriptions with the Maxwell equations leads to an eigenvalue problem for each layer and finally to the propagating and evanescent modes inside the layers. By applying the appropriate boundary conditions all layers are coupled and the resulting reflected and transmitted plane waves at the top an bottom side of the mask (see Figure 1) are computed. Because of the periodic boundary conditions in lateral direction, the mask is always regarded as a periodic structure. Isolated features can be simulated by using a mask period which is large enough for the respective problem. Two dimensional (e.g. lines) and three dimensional (e.g. contact holes) computations are possible with the method.

The computation time is proportional to M3. M is the over-all number of modes taken into account for the computation. Except of very special mask geometries (e.g. large EUV mask with defective and deformed multilayer) and assuming that up-to-date computer systems have usually one or more GB of working memory the memory requirement of the method is small and no limiting factor.

The method is qualified for lithography simulations and in particular for EUV lithography simulations because of the following reasons:

  • Good convergence for all mask types investigated so far
  • Mask geometry usually corresponds to the required slicing
  • Only stationary fields are of interest, no simulations in the time domain required, the waveguide method is operating in the frequency domain → fast simulations
  • Computation time depends on the used number M of modes, M usually low in case of typical simulation geometries → fast simulations
  • In the EUV case: Refractive indices of all materials close to 1, reduced number of modes required for the field and material description → fast EUV simulations
  • In the two dimensional case fast simulations even for a high number M of used modes

In the optical case (e.g. at 193 nm), the convergence is a basic problem of all electromagnetic field solvers based on modal methods like waveguide because of the significant refractive index differences between the materials. An implemented mathematical optimization reduces this problem effectively with the result of a very good convergence in this case. In contrast to that the EUV case is uncritical.

Defective EUV masks

A special feature of both simulation methods implemented at Fraunhofer IISB is the capability to deal with EUV masks with a defective multilayer. Figure 2 shows the simulated cross section of a multilayer with a Gaussian shaped defect onto the substrate.

Figure 2: Simulated cross section of a defective EUV multilayer


The complete topography of the defect and of the deformed multilayer can be characterized by many parameters. In the example of Figure 2 a Gaussian shaped multilayer defect causes a deformation of the whole multilayer structure. Other shapes can be added easily to the simulation model if corresponding measurement data are available. If other data concerning the layer deformation are available new compression models can be added easily to the simulation system.

With the waveguide method the whole multilayer is simulated rigorously in the same way as described above. In case of three dimensional simulations a defective multilayer geometry can have a significant impact on the computation time. In case of two dimensional simulations the impact on the computation time is small due to the fact that two dimensional waveguide simulations are very fast in general.

With the FDTD method the mask is split up into the absorber part which is computed rigorously and in the multilayer part which is computed analytically based on a transfer matrix method.

Decomposition technique

If larger mask areas have to be simulated the computation time of a three dimensional simulation can become too long. A simplified three dimensional computation based on a decomposition technique can be applied. This technique replaces a full three dimensional simulation by a superposition of several two dimensional and one dimensional simulations. The result is a significantly reduced computation time but compared to a full three dimensional simulation an error must be accepted.

The basic idea of the decomposition technique is to separate diffraction effects from mask edges along x- and y-direction. A full three dimensional electromagnetic field simulation is split into several two dimensional and one dimensional simulations. The decomposition technique is implemented for the FDTD and for the waveguide method.
Figure 3: Topview on a three dimensional EUV mask, decomposition into two dimensional cuts (cut 1-6) with constant dielectric properties in x- and y-direction, respectively

According to the example of Figure 3 three parts in x-direction (denoted by part 1, part 2, part 3) and 3 parts in y-direction (denoted by part 4, part 5, part 6) can be found. All two dimensional cuts representing the mask parts 1 to 6 (denoted by cut 1 to cut 6) are computed rigorously under the assumption that the dielectric properties of the mask vary only in x- and y-direction respectively. The resulting field (FDTD simulations) or spectrum (waveguide simulations) of the three dimensional mask is obtained by a composition of the complex fields (FDTD simulations ) or spectra (waveguide simulations) of the x- and y-configurations and the mask reflectivities. For a better performance the spectrum composition of a waveguide computation is performed in the frequency domain.

Simulation examples

(If not mentioned otherwise all dimensions are in mask scale)

Figure 4: Waveguide simulations: Mask spectrum simulation time depending on the number of used diffraction orders for 3D simulations (left picture) and 2D simulations (right picture) on a 1.9 GHz laptop. The indicated number N of orders is not the above mentioned over-all number of orders. The over-all number M of orders is (2N + 1) 2 (3D simulation, left picture) and 2N+1 (2D simulation, right picture). 3D simulations: N = 17 corresponds to a simulation area of about 1.2 µm x 1.2 µm and a computation time of about 4600 s, 2D simulations: N = 500 corresponds to a simulation area of about 32.2 µm and a computation time of about 314 s. The computation times e.g. for a 90 nm contact and a 90 nm line (both wafer scale) are 337 s and 0.1 s respectively. In general the size b of the simulation area in x-direction (2D simulations) and in x- and y-direction (3D simulations) is b = N * λ / 3 for optical masks (e.g. at 193 nm) and b = N * λ * 2 for EUV masks (at 13.5 nm). λ is the wavelength.


Figure 5: Waveguide simulations: Simulated mask near fields of two large EUV mask areas with several post absorbers onto a typical EUV multilayer. The simulation area of the left picture is 100 λ * 100 λ and of the right picture 200 λ * 200 λ. λ is the wavelength. The application of the decomposition technique in combination with waveguide leads to a mask near field simulation time of only 20 s for the left hand picture and 220 s for the right hand picture on a 1.9 GHz laptop. Another application of the decomposition technique is the simulation of masks with standard sizes to get quickly a first result. For instance the simulation time of an optical mask with a 90 nm contact (pitch 180 nm, wafer scale) is only 0.2 s compared to 337 s for the same real three dimensional computation. It is important to mention that the application of the decomposition technique is only an approximation of a real three dimensional computation. However, at the moment the decomposition technique is the only possibility to simulate large mask areas rigorously.
Figure 6: Waveguide simulations: Computed mask near fields of a post absorber onto a typical EUV multilayer with a multilayer defect at different positions with respect to the absorber. The rectangular dark parts in the pictures result from the rectangular post absorber, the roundish dark areas from the multilayer defect. Depending on the defective multilayer structure the simulation time can be significantly longer compared to standard masks. For the shown examples the simulation time was about 9000 s.
Figure 7: FDTD simulation: Time resolved simulation of light diffraction at an EUV mask with a line absorber, upper line from left to right: Downward propagating incident light (the blue part shows the area of the line absorber), lower line from left to right: After reflection at the multilayer upward propagating reflected light (the blue part shows the area of the line absorber). The typical standing waves can be seen.