3D Simulations of Electromagnetic Fields in Nanostructures using the Time-Harmonic Finite-Element Method

3D Simulations of Electromagnetic Fields in Nanostructures using the Time-Harmonic Finite-Element Method


Rigorous computer simulations of propagating electromagnetic fields have become an important tool for optical metrology and optics design of nanostructured components. As has been shown in previous benchmarks some of the presently used methods suffer from low convergence rates and/or low accuracy of the results and exhibit very long computation times [1, 2] which makes application to extended 2D layout patterns impractical. We address 3D simulation tasks by using a finite-element solver which has been shown to be superior to competing methods by several orders of magnitude in accuracy and computational time for typical microlithography simulations [2]. We report on the current status of the solver, incorporating higher order edge elements, adaptive refinement methods, and fast solution algorithms. Further, we investigate the performance of the solver in the 3D simulation project of light diffraction off an alternating phase-shift contact-hole mask.

3D electromagnetic simulations, microlithography, finite-element methods, FEM

Corresponding author: S. Burger
URL: http://www.zib.de/Numerik/NanoOptics/
Email: burger@zib.de Copyright 2007 Society of Photo-Optical Instrumentation Engineers.
This paper will be published in Proc. SPIE Vol. 6617 (2007), (Modeling Aspects in Optical Metrology, H. Bosse, Ed.) and is made available as an electronic preprint with permission of SPIE. One print or electronic copy may be made for personal use only. Systematic or multiple reproduction, distribution to multiple locations via electronic or other means, duplication of any material in this paper for a fee or for commercial purposes, or modification of the content of the paper are prohibited.

1 Introduction

With the advances of micro- and nanotechnology simulation tools for rigorous solutions of Maxwell’s equations have become an important tool in research and development. Designing, e.g., a nanooptical component or a metrology tool is usually assisted by computer simulations, and simulations are used in most scientific works on nanooptical research to support theoretical and experimental findings. It exists a variety of different methods for solving Maxwell’s equations, and generally also a variety of different numerical implementations of each method. Prominent examples of different methods are the finite-element method (FEM), the finite difference time domain method (FDTD), the boundary element method (BEM), and rigorously coupled wave analysis (RCWA).

We address 3D simulation tasks occuring in microlithography by using the frequency-domain finite-element solver JCMsuite. This solver has been successfully applied to a wide range of 3D electromagnetic field computations including microlithography [1, 2, 3], left-handed metamaterials in the optical regime [4, 5], photonic crystals [6], and nearfield-microscopy [7]. The solver has also been used for pattern reconstruction in EUV scatterometry [8, 9], and it has been benchmarked to other methods in typical DUV lithography [1, 2] and other [10, 6] projects.

In this paper we report on the current status of the finite-element solver JCMsuite, and we present simulations of light transition through periodic arrays of alternating phase-shift contact-holes which have to be printed on wafers with the same reticle with high overlapping process window. 3D effects cause imbalances of intensities of the - and -holes, and of different behavior in defocus if the etch depth is not adjusted. The evaluation of the impact of these effects has been chosen as a test case for JCMsuite.

2 Background

The finite-element package JCMsuite allows to simulate a variety of electromagnetic problems. It incorporates a scattering solver (JCMharmony), a propagating mode solver (JCMmode) and a resonance solver (JCMresonance). The scattering, eigenmode and resonance problems can be formulated on 1D, 2D and 3D computational domains. Admissible geometries can consist of periodic or isolated patterns, or a mixture of both. Further, solvers for problems posed on cylindrically symmetric geometries are implemented.


0d \psfrag1d \psfragpx \psfragpy

Figure 1: Schematics of the 3D test structure: cross-section in a --plane. 0-degree (180-degree) rectangular phase-shift holes are depicted in grey (light grey). The elementary vectors, , of the periodic pattern are indicated by arrows.

p_x \psfrag0d \psfrag1d \psfragd3 \psfragd2 \psfragd1 \psfragdetch \psfragsubstratesubstrate \psfragsuperspacesuperspace (air)

Figure 2: Schematics of the 3D test structure: cross-section in a --plane. The material stack on a Quartz (SiO) substrate consists of a molydenum silicide and two layers containing chromium. 0-degree (180-degree) rectangular phase-shift air-filled holes are indicated.

cdx0 \psfragcdx1 \psfragcdy0 \psfragcdy1 \psfragcomp_domaincomputational domain

Figure 3: Schematics of the 3D test structure: cross-section in a --plane. The computational domain is indicated by a polygon.

In this paper we concentrate on light-scattering off a 3D pattern (photomask) which is periodic in the and directions and is enclosed by homogeneous substrate (at ) and superstrate (at ) which are infinite in the , resp. direction (see Figures 1 and 2). Light propagation in the investigated system is governed by Maxwell’s equations where vanishing densities of free charges and currents are assumed [11]. The dielectric coefficient and the permeability of the considered photomasks are periodic and complex, , . Here is any elementary vector of the periodic lattice. For given primitive lattice vectors and an elementary cell can be defined as . As can be seen from the computational domain in Figure 3 elementary cells of the same volume can be defined also by more complicated polygons which are adapted to the actual geometry. 1

A time-harmonic ansatz with frequency and magnetic field leads to the following equations for :

  • The wave equation and the divergence condition for the magnetic field:

  • Transparent boundary conditions at the boundaries to the substrate (at ) and superstrate (at ), , where is the incident magnetic field (plane wave in this case), and is the normal vector on :


    The operator (Dirichlet-to-Neumann) is realized with an adaptive PML method [12, 13]. This is a generalized formulation of Sommerfeld’s radiation condition.

  • Periodic boundary conditions for the transverse boundaries, , governed by Bloch’s theorem [14]:


    where the Bloch wavevector is defined by the incoming plane wave .

Similar equations are found for the electric field ; these are treated accordingly. The finite-element method solves Eqs. (1) – (4) in their weak form, i.e., in an integral representation.

The finite-element methods consists of the following steps:

  • The computational domain is discretized with simple geometrical patches, JCMsuite uses linear (1D), triangular (2D) and tetrahedral or prismatoidal (3D) patches. The use of prismatoidal patches is well suited for layered geometries, as in photomask simulations.

  • The function spaces in the integral representation of Maxwell’s equations are discretized using Nedelec’s edge elements, which are vectorial functions of polynomial order defined on the simple geometrical patches [15]. In the current implementation, JCMsuite uses polynomials of first () to ninth () order. In a nutshell, FEM can be explained as expanding the field corresponding to the exact solution of Equation (1) in the basis given by these elements.

  • This expansion leads to a large sparse matrix equation (algebraic problem). To solve the algebraic problem on a standard workstation linear algebra decomposition techniques (LU-factorization, e.g., package PARDISO [16], which was used in the simulations of Chapter 3) are used. In cases with either large computational domains or high accuracy demands, also domain decomposition methods [17] are used and allow to handle problems with very large numbers of unknowns.

For details on the weak formulation, the choice of Bloch-periodic functional spaces, the FEM discretization, and the implementation of the adaptive PML method in JCMsuite we refer to previous works [12, 13]. In future implementations performance will further be increased by using curvilinear elements, general domain-decomposition techniques and -adaptive methods.

3 3D simulations of a contact-hole photomask

3.1 Parameter settings and simulation flow

We apply the 3D light-scattering solver of the programme package JCMsuite to the test case of a 2D periodic pattern of an alternating phase-shift contact-hole photomask. Figures 12 and 3 show cross-sections through the geometrical layout of the mask and define the geometrical parameters. The geometrical parameters investigated in this study are given in Table 1. The computational domain is a unit-cell of the periodic pattern; it is indicated in Figure 3. Its total volume is , where . Please note that a simple rectangular box-like computational domain would have twice the volume, . Within the numerical error the results are independent of the choice of the unit-cell chosen as computational domain.


LayoutGeometry layout \psfragTriangulationMesh parameters \psfragMaterialsMaterial parameters \psfragMeshMesh \psfragSourcesSources parameters \psfragProjectProject parameters \psfragResultsResults \psfragProjectionProjection \psfragInterfaceMATLAB Interface: Automatic input generation, parameter scans, data analysis

Figure 4: Schematics of the simulation flow of the FEM solver JCMsuite.
parameter set 1 set 2 set 3 set 4 set 5 set 6
[nm] 576 640 720 800 1040 1200
[nm] 1120 1120 1120 1120 1120 1120
[nm] 416 452 496 552 496 520
[nm] 876 800 752 696 976 980
3D bias [nm] (32, 48, 64, 80)
[nm] 20
[nm] 40
[nm] 68
[nm] (156, 160, 164, 168, 172, 176, 180)
193.0 nm
Table 1: Parameter settings for the 3D simulations. The relative permittivity, , is the square of the complex index of refraction, . The simulation wavelength was nm while the wavelength in the real system is nm. The dispersion of is taken into account at nm. All combinations of 3D bias and etch depth, , have been investigated for the six parameter sets in order to find for each pattern the best value and the best overall 3D bias and etch depth.

(a)(a) \psfrag(b)(b) \psfrag(c)(c)

Figure 5: Cross-sections through prismatoidal discretizations of several computational domains with different geometry parameters, corresponding to parameter sets 1 (a), 4 (b), 6 (c)(see Table 1). Please note that (a), (b), and (c) have different length scales.

Figure 4 shows the simulation flow of the finite-element software:

  • The geometry of the computational domain is described in a polygonal format of the cross-section, corresponding to Fig. 1, and includes a height profile (see Fig. 2) with material attributions and statements about the computational domain boundaries (periodic and transparent boundaries in this case). The translational vectors of the periodic pattern () are identified automatically from the layout, optimized settings of the perfectly matched layers (PML) are found automatically with an adaptive method  (adaptive PML, aPML). Further, parameters specifying a maximum patch size of the finite-element mesh and further meshing properties can be set. From these geometry parameters and mesh parameters the prismatoidal mesh is generated automatically. Figure 5 shows -cross-sections through the prismatoidal meshes for different geometrical parameter sets.

  • Material parameters (complex permittivity and permeability tensors) can be specified as piecewise constant functions and/or (using dynamically loaded libraries) as functions of arbitrary spatial dependence. In the presented example, the piecewise constant, isotropic settings given in Table 1 are used. The relative permeability is for all used materials.

  • Light sources can be defined as predefined functions (plane waves, Gaussian beams, point sources) or as arbitrary functions using dynamically loaded libraries. JCMsuite allows to generate solutions to several independent source terms in parallel, efficiently re-using the inverted system matrix. In the presented example, two plane waves with normal incidence (wavevector ), and with orthogonal polarizations are used.

  • The main project definitions in this case are the accuracy settings (mesh refinement, PML refinement, finite-element degree) and the definitions of postprocesses. Upon execution the JCMsolve computes the full field distribution over the entire computational domain. Through internal or external post-processes, the quantities of interest can be derived from the field. E.g., the complex amplitudes of propagating modes are attained by Fourier transformation, an aerial image is calculated [3], or the field distribution is exported to graphics format for visualisation and analysis.

  • Interfaces to scripting languages like Python or MATLAB can be used for performing automatic parameter scans and data analysis. In the presented example, the MATLAB interface for the generation of the input files and for the execution of a loop over the geometrical parameters given in Table 1 consists of well below 100 lines of MATLAB code.

3.2 Error analysis


(a)(a) \psfrag(b)(b)

Figure 6: Convergence of the intensities of the first (a) and zero (b) diffraction orders with increasing number of unknowns of the FEM solution. Results from computations with finite element degree are displayed.

(a)(a) \psfrag(b)(b)

Figure 7: (a) Convergence of the total transmission amplitude to the 5 x 5 lowest diffraction orders. (b) Computation times in dependence of the number of unknowns of the corresponding matrix equation.

We have investigated the accuracy of the near field results by simulating near fields using finite-elements of different polynomial degrees (-refinement) and using meshes with different numbers of mesh elements (-refinement) for one fixed setting of geometrical and other physical parameters. From each near fields we compute the diffraction pattern using a post-process. Figure 6 shows the convergence of the absolute values of the intensities in the zero and first diffraction orders. Figure 7(a) shows the convergence of the relative error of the intensity in the 5 x 5 lowest diffraction orders. The relative error of the intensity, , is defined as , where the quasi-exact intensity, , is deduced from a solution with highest finite-element degree and finest mesh. As can be seen from Figure 7, relative errors of below one percent for the prominent diffraction orders can be reached using a high polynomial degree of the finite-elements () and relatively moderate numbers of unknowns, corresponding to coarse spatial discretizations.

3.3 Computational costs

Figure 7(b) shows the computation time in minutes in dependence on the number of degrees of freedom in the simulation. Different symbols in the plot correspond to different finite-element polynomial degrees. The computations have been performed on a standard workstation with 16 (2 x 8) processors and extended RAM (AMD Opteron, 64GB RAM). Please note that for the results in Figure 7(b) we have used only four of the 16 processors, while for the scans over geometrical parameters we have used twelve out of the 16 processors, giving a speed-up factor of approximately three in the computation times. As can be seen in Figure 6 and Figure 7, results corresponding to an accuracy of the prominent diffraction orders of better 1% can be attained in computation times of few minutes when using finite elements of order and above. Memory requirements were roughly 5…15GB RAM for and finite-element degree (the setting used for the scan over the geometrical parameters). Please note that in general it does depend on the physical problem wheather a refinement of the finite element degree , , or a refinement of the mesh width , , leads to a better convergence of the results. For the given rather large (on the wavelength-scale) structures on the photomask, rather high degrees seem to lead to the best results, while for problems with smaller features finer meshes are favourable. In a future implementation of the programme package this choice will be automatically and locally detected by -adaptive methods.


(a)(a) \psfrag(b)(b)

Figure 8: Intensity of the first (+1,0) diffraction order for the different investigated parameter sets. Each parameter set consists of four groups (3D bias=8,12,16,20 nm) of seven values of etch depths ( nm).

NN [  ] \psfragapprox. total time [min]total time [ min ]

Figure 9: Intensity of the zero (0,0) diffraction order for the different investigated parameter sets. Each parameter set consists of four groups (3D bias=8,12,16,20 nm) of seven values of etch depths ( nm).

3.4 Geometry dependence of the diffraction spectrum

We have performed a scan over the geometrical parameters defined in Table 1. The six parameter sets consist each of 28 subsets of the various combinations of 3D bias and etch depth. Figure 8 and 9 show how the first (+1,0) order and zero order diffraction intensities vary over the scan of geometrical parameters. As expected for phase-shift masks, the intensity in the first diffraction orders (prominent order) is about two orders of magnitude higher than the intensity in the zero diffraction order. For otherwise fixed geometrical parameters, the variation of the prominent diffraction order with the etch depth in the investigated regime is about 2-4%. The variation of the prominent diffraction order with the 3D bias is about 10-15%. The variation with a different set of CD and pitch parameters is much larger.

This shows that for mask design using rigorous simulations, especially for the design of parameters such as the etch depth, it is necessary to use simulation tools which produce results of low numerical error.


(a)(a) \psfrag(b)(b)

Figure 10: (a) CD as a function of defocus for an almost balanced case. (b) CD as a function of defocus for an imbalanced case.

3.5 Imaging application

For the conditions of Table 1 the farfield coefficients were determined and used as input for an Qimonda inhouse imaging simulator. The mask bias for each pitch was pre-optimized by Kirchhoff (2D) mask simulations in resist. For evaluation of the arial image output for each pitch and 3D bias the intensity threshold for target CD was determined at best focus as average of 0- and 180-holes. Figure 10(b) shows the typical phase balancing problem that occurs in case of non-optimum etch depth in the mask: the CD of 0- and 180-holes behaves different at defocus, whereas for the optimum depth the behavior is symmetric (Figure 10(a)).

The CD difference of 0- and 180-holes over defocus shows a slope that is characteristic for the phase error, i.e., the deviation from optimum etch depth where the defocus behavior is flat (see Figure 11(a)). Furtheron, the 180-hole prints smaller than a 0-hole of same size due to 3D mask effects. This is compensated by a 3D mask bias, see Figures 11(b) and 12(a), which is in the range of 8 – 20 nm, depending on pattern and mask conditions (side wall angles, material). The slope is rather independent of the 3D bias, as shown in Figure 12(b). The optimum 3D bias, where the CD difference is zero, is 8.5 nm at an etch depth of 168 nm, as obtained by a linear fit of the CD difference vs. 3D bias. As can be also seen in Figure 12(b), the slope is almost zero at this etch depth.


(a)(a) \psfrag(b)(b)

Figure 11: (a) CD difference () as a function of defocus and etch depth for one fixed 3D bias. (b) CD difference () as a function of defocus and 3D bias at optimized etch depth.

A determination of the optimum etch depths and optimum 3D mask biases through pitch reveals significant differences. While the 3D mask bias can be corrected individually for such special patterns by software algorithms it is challenging to correct the 3D mask effect in real layouts containing random patterns. An even more critical issue is to find an optimum etch depth through pitch. Here an analysis of overlapping process windows including etch depth deviations is the next step to find an optimum etch depth through pitch. Then conclusions can be drawn with respect to lithographic applications for nodes below 40 nm half pitch.


(a)(a) \psfrag(b)(b)

Figure 12: (a) CD difference () as a function of 3D bias at optimized etch depth (168 nm). (b) Slope of CD difference, , as a function of etch depth and 3D bias.

4 Conclusions

We have performed rigorous 3D FEM simulations of light transition through alternating phase-shift contact-hole photomasks using the finite-element programme package JCMsuite. We have investigated the convergence behavior of the solutions and we have shown that we achieve results at high numerical accuracy. Our results show that rigorous 3D mask simulations can well be handled at high accuracy and relatively low computational cost.

The investigation of the diffraction spectrum of an alternating phase-shift contact-hole mask has shown that for mask design and optimization projects it is necessary to use rigorous tools with low discretization errors, typically below 1% for the prominent diffraction orders, because the variance of the diffraction orders can be of the order of few percent for typical design parameter regimes. We have used the simulations to determine optimum etch depths and 3D mask biases for the various investigated parameter sets.

This paper was supported by the EFRE fund of the European Community and by funding of the State Saxony of the Federal Republic of Germany (project number 10834). The authors are responsible for the content of the paper.


  1. In this case a computational domain constructed corresponding to would in general lead to intersections of the domain boundary and geometrical features at sharp angles. This could result in low quality triangulations.


  1. S. Burger, R. Köhle, L. Zschiedrich, W. Gao, F. Schmidt, R. März, and C. Nölscher, “Benchmark of FEM, waveguide and FDTD algorithms for rigorous mask simulation,” in Photomask Technology, J. T. Weed and P. M. Martin, eds., 5992, pp. 378–389, Proc. SPIE, 2005.
  2. S. Burger, R. Köhle, L. Zschiedrich, H. Nguyen, F. Schmidt, R. März, and C. Nölscher, “Rigorous simulation of 3d masks,” in Photomask Technology, P. M. Martin and R. J. Naber, eds., 6349, p. 63494Z, Proc. SPIE, 2006.
  3. R. Köhle, “Rigorous simulation study of mask gratings at conical illumination,” Proc. SPIE 6607, pp. 6607–106, 2007.
  4. S. Linden, C. Enkrich, G. Dolling, M. W. Klein, J. Zhou, T. Koschny, C. M. Soukoulis, S. Burger, F. Schmidt, and M. Wegener, “Photonic metamaterials: Magnetism at optical frequencies,” IEEE Journal of Selected Topics in Quantum Electronics 12, pp. 1097–1105, 2006.
  5. G. Dolling, M. Wegener, A. Schädle, S. Burger, and S. Linden, “Observation of magnetization waves in negative-index photonic metamaterials,” Appl. Phys. Lett. 89, p. 231118, 2006.
  6. S. Burger, R. Klose, A. Schädle, F. Schmidt, and L. Zschiedrich, “Adaptive FEM solver for the computation of electromagnetic eigenmodes in 3d photonic crystal structures,” in Scientific Computing in Electrical Engineering, A. M. Anile, G. Ali, and G. Mascali, eds., pp. 169–175, Springer Verlag, 2006.
  7. T. Kalkbrenner, U. Hakanson, A. Schädle, S. Burger, C. Henkel, and V. Sandoghdar, “Optical microscopy using the spectral modifications of a nano-antenna,” Phys. Rev. Lett. 95, p. 200801, 2005.
  8. J. Pomplun, S. Burger, F. Schmidt, L. W. Zschiedrich, F. Scholze, and U. Dersch, “Rigorous FEM-simulation of EUV-masks: Influence of shape and material parameters,” in Photomask Technology, 6349-128, Proc. SPIE, 2006.
  9. Y. Tezuka, J. Cullins, Y. Tanaka, T. Hashimoto, I. Nishiyama, and T. Shoki, “EUV exposure experiment using programmed multilayer defects for refining printability simulation,” 6517, Proc. SPIE, 2007.
  10. R. Holzlöhner, S. Burger, P. J. Roberts, and J. Pomplun, “Efficient optimization of hollow-core photonic crystal fiber design using the finite-element method,” J. Europ. Opt. Soc: Rap. Comm. 1, p. 06011, 2006.
  11. A. K. Wong, Optical Imaging in Projection Microlithography, SPIE Press, Bellingham, 2005.
  12. L. Zschiedrich, R. Klose, A. Schädle, and F. Schmidt, “A new finite element realization of the Perfectly Matched Layer Method for Helmholtz scattering problems on polygonal domains in 2D,” J. Comput. Appl. Math. 188, pp. 12–32, 2006.
  13. L. Zschiedrich, S. Burger, B. Kettner, and F. Schmidt, “Advanced finite element method for nano-resonators,” 6115, p. 41, Proc. SPIE, 2006.
  14. K. Sakoda, Optical Properties of Photonic Crystals, Springer-Verlag, Berlin, 2001.
  15. P. Monk, Finite Element Methods for Maxwell’s Equations, Claredon Press, Oxford, 2003.
  16. O. Schenk et al., “Parallel sparse direct linear solver PARDISO.” Department of Computer Science, Universität Basel.
  17. L. Zschiedrich, S. Burger, A. Schädle, and F. Schmidt, “Domain decomposition method for electromagnetic scattering problems,” in Proceedings of the 5th International Conference on Numerical Simulation of Optoelectronic devices, pp. 55–56, 2005.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description