3D Simulations of Electromagnetic Fields in Nanostructures using the TimeHarmonic FiniteElement Method
Abstract
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 finiteelement 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 phaseshift contacthole mask.
Corresponding author: S. Burger
URL: http://www.zib.de/Numerik/NanoOptics/
Email: burger@zib.de
Copyright 2007 Society of PhotoOptical 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 finiteelement 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 frequencydomain finiteelement solver JCMsuite. This solver has been successfully applied to a wide range of 3D electromagnetic field computations including microlithography [1, 2, 3], lefthanded metamaterials in the optical regime [4, 5], photonic crystals [6], and nearfieldmicroscopy [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 finiteelement solver JCMsuite, and we present simulations of light transition through periodic arrays of alternating phaseshift contactholes 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 finiteelement 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.
In this paper we concentrate on
lightscattering 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.
A timeharmonic ansatz with frequency and magnetic field leads to the following equations for :

The wave equation and the divergence condition for the magnetic field:
(1) (2) 
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 :
(3) The operator (DirichlettoNeumann) 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]:
(4) where the Bloch wavevector is defined by the incoming plane wave .
Similar equations are found for the electric field ; these are treated accordingly. The finiteelement method solves Eqs. (1) – (4) in their weak form, i.e., in an integral representation.
The finiteelement 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 (LUfactorization, 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 Blochperiodic 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 domaindecomposition techniques and adaptive methods.
3 3D simulations of a contacthole photomask
3.1 Parameter settings and simulation flow
We apply the 3D lightscattering solver of the programme package JCMsuite to the test case of a 2D periodic pattern of an alternating phaseshift contacthole photomask. Figures 1, 2 and 3 show crosssections 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 unitcell of the periodic pattern; it is indicated in Figure 3. Its total volume is , where . Please note that a simple rectangular boxlike computational domain would have twice the volume, . Within the numerical error the results are independent of the choice of the unitcell chosen as computational domain.
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 
[nm]  
[nm]  
3D bias [nm]  (32, 48, 64, 80)  
[nm]  20  
[nm]  40  
[nm]  68  
[nm]  (156, 160, 164, 168, 172, 176, 180)  
193.0 nm  
1.0  
Figure 4 shows the simulation flow of the finiteelement software:

The geometry of the computational domain is described in a polygonal format of the crosssection, 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 finiteelement mesh and further meshing properties can be set. From these geometry parameters and mesh parameters the prismatoidal mesh is generated automatically. Figure 5 shows crosssections 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 reusing 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, finiteelement degree) and the definitions of postprocesses. Upon execution the JCMsolve computes the full field distribution over the entire computational domain. Through internal or external postprocesses, 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
We have investigated the accuracy of the near field results by simulating near fields using finiteelements 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 postprocess. 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 quasiexact intensity, , is deduced from a solution with highest finiteelement 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 finiteelements () 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 finiteelement 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 speedup 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 finiteelement 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 wavelengthscale) 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.
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 phaseshift 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 24%. The variation of the prominent diffraction order with the 3D bias is about 1015%. 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.
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 preoptimized 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 180holes. Figure 10(b) shows the typical phase balancing problem that occurs in case of nonoptimum etch depth in the mask: the CD of 0 and 180holes behaves different at defocus, whereas for the optimum depth the behavior is symmetric (Figure 10(a)).
The CD difference of 0 and 180holes 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 180hole prints smaller than a 0hole 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 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.
4 Conclusions
We have performed rigorous 3D FEM simulations of light transition through alternating phaseshift contacthole photomasks using the finiteelement 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 phaseshift contacthole 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.
Acknowledgements.
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.Footnotes
 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.
References
 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.
 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.
 R. Köhle, “Rigorous simulation study of mask gratings at conical illumination,” Proc. SPIE 6607, pp. 6607–106, 2007.
 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.
 G. Dolling, M. Wegener, A. Schädle, S. Burger, and S. Linden, “Observation of magnetization waves in negativeindex photonic metamaterials,” Appl. Phys. Lett. 89, p. 231118, 2006.
 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.
 T. Kalkbrenner, U. Hakanson, A. Schädle, S. Burger, C. Henkel, and V. Sandoghdar, “Optical microscopy using the spectral modifications of a nanoantenna,” Phys. Rev. Lett. 95, p. 200801, 2005.
 J. Pomplun, S. Burger, F. Schmidt, L. W. Zschiedrich, F. Scholze, and U. Dersch, “Rigorous FEMsimulation of EUVmasks: Influence of shape and material parameters,” in Photomask Technology, 6349128, Proc. SPIE, 2006.
 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.
 R. Holzlöhner, S. Burger, P. J. Roberts, and J. Pomplun, “Efficient optimization of hollowcore photonic crystal fiber design using the finiteelement method,” J. Europ. Opt. Soc: Rap. Comm. 1, p. 06011, 2006.
 A. K. Wong, Optical Imaging in Projection Microlithography, SPIE Press, Bellingham, 2005.
 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.
 L. Zschiedrich, S. Burger, B. Kettner, and F. Schmidt, “Advanced finite element method for nanoresonators,” 6115, p. 41, Proc. SPIE, 2006.
 K. Sakoda, Optical Properties of Photonic Crystals, SpringerVerlag, Berlin, 2001.
 P. Monk, Finite Element Methods for Maxwell’s Equations, Claredon Press, Oxford, 2003.
 O. Schenk et al., “Parallel sparse direct linear solver PARDISO.” Department of Computer Science, Universität Basel.
 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.