Open Access
Add to BibSonomyAbstract
Early formulations of the RCWA yield, implicated by the erroneous application of factorization rules to discrete Fourier transformations, poor convergence in certain cases. An explanation for this finding and an approach to overcome the problem for crossed gratings was first given by Li [J. Opt. Soc. Am. A 13, 1870 (1996) and 14, 2758 (1997)]. A further improvement was achieved by Schuster et al. [J. Opt. Soc. Am. A 24, 2880 (2007)], using a structure dependent normal vector (NV) field. While it is trivial to create those NV fields for simple geometrical shapes, to our knowledge an appropriate algorithm for arbitrary shapes does not exist, yet. In this work we present such an algorithm.
©2008 Optical Society of America
1. Introduction
For the numerical simulation of grating diffraction rigorous coupled-wave analysis (RCWA) is a commonly used method. After its publication in 1981 by Moharam and Gaylord [1] it turned out that its results for metallic gratings in TM polarization suffered from slow convergence [2] due to numerical problems. While empirical solutions came from Lalanne [3] as well as from Guizal and Granet [4], Li finally gave a theoretical explanation, introducing his factorization rules [5]. Since then different approaches were pursued to apply these rules to crossed gratings in order to improve convergence [6, 7]. One of the most recent improvements for crossed gratings was presented by Schuster et al. [8]. They adopted the idea of the Fast Fourier Factorization by Nevière and Popov [9, 10] and developed a normal vector (NV) method which is an optimized solution for the correct application of Li’s factorization rules to crossed gratings. This method, however, requires a structure dependent NV field. The requirements for this vector field will be explained in the next Section (for a detailed explanation see [8]). Afterward, Section 3 and 4 will explain how the algorithm works. Results of the algorithm will be shown in Section 5.
2. Background
The RCWA is based on a time-harmonic version of Maxwell’s Equations to solve the diffraction problem:
To solve these equations the RCWA uses pseudo Fourier expansions as an ansatz for E, H and ε to transform this system of PDEs to a system of ODEs, which, in turn, can be reduced to an eigenvalue problem. Important for the following considerations is the calculation of the product ε0ε E. This product, the dielectric displacement
transforms into a convolution in Fourier space. At material boundaries the normal component of Eq. (3) – in contrast to the tangential component – is a product of complementary discontinuous functions ε and E ⊥, as is generally known from electrodynamics. In truncated Fourier space, as it is used in numerical calculations, this product cannot be calculated by the usual convolution
as this leads to unnecessary oscillations in the dielectric displacement D ⊥ as shown in Fig. 4 of [5]. In Eq. (4) single brackets refer to vertical vectors containing all Fourier coefficients of the specified function and 〚ε〛 denotes the Toeplitz matrix of the Fourier coefficients of ε with the entries εjk=εj-k.
While Lalanne gave strong empirical evidence [3] that another rule must be used for the normal component of the electrical field, Li later gave a theoretical explanation [5] why the so-called inverse rule
must be used for complementary discontinuous functions instead.
Consequently the normal and tangential components must be treated differently. While Li used a method of zigzag lines [7] and got an improved convergence, an even further improvement in convergence could be achieved by the normal vector method of Schuster et al. [8].
This method basically projects the electrical field onto normal vectors to get the tangential and normal components of the electrical field. Until now, however, there was no way to calculate those normal vector fields for completely arbitrary structures. A solution for such an algorithm will be given in the following.
3. Automatic generation of normal vector fields
The algorithm is based on the fact that normal vectors are only defined on the boundaries. Except for their orientation (inwards or outwards) vectors are unambiguously defined here. Since the products N2x, N2y and NxNy are transformed into Fourier space, as can be seen in Eq. (8) of [8], the criterion for the vectors in the remaining space is smooth transitions between NVs. An interpolation fulfills this requirement. Generating the vector field is consequently done in two steps:
1. Calculating the normal vectors on the boundaries, called boundary vectors
2. Interpolating the vectors in the remaining space
These two steps will be explained further in the following two subsections.
3.1. Calculating the boundary vectors
A bitmap serves as a description for the crossed grating structure, where different materials are indicated by different color indices. Determining the boundary vectors is done by calculating the gradient of this bitmap. The pixelation yields a staircase-like modulation of the boundary vectors (Fig. 1a). A blurring of the bitmap by a Gaussian filter smooths the pixelation and hence avoids the modulation (Fig. 1b). However, the filter also smooths the distinct jumps in the bitmap and hence results in a wide band of normal vectors, which is undesired. The gradient of the original image Fig. 1a is only ≠0 in a narrow band around the boundary and thus can be used as a mask for Fig. 1b. The result of this procedure is a narrow band of proper boundary vectors as can be seen in Fig. 1c. After the boundary vectors are calculated the orientation of the vectors can be given a preferred orientation as it is shown in Fig. 1d. In this way a degree of freedom is introduced which causes differences in both the appearance of the NV fields and the end results of the RCWA. The influence of a preferred orientation is investigated in Section 5.
Fig. 1. (a) Staircase effect of the gradient (b) Gradient after blurring the bitmap (c) Result of combining (a) and (b) (d) Optional change of the preferred direction
3.2. Calculating the remaining vectors
Interpolation is done using inverse distance weighting [11, 12]. This is an interpolation method suited for irregularly spaced nodes which can be implemented to work very fast as will be shown in the next Section. Using this interpolation method, a vector N* at position x can be calculated as
and be normalized as
.
N i denotes the boundary vectors, x i their positions and nb their count. Eq. (6) must be evaluated for every x except for those where boundary vectors already exist.
4. Accelerating the algorithm
Using this implementation of the algorithm one realizes that the time complexity is 𝒪(R3), where R×R is the resolution. With a resolution of e.g. 2048×2048 this leads to very long computation durations. To gain an advantage for the NV method against other methods the complexity can be reduced to 𝒪(R2) using a progressive refinement algorithm. This algorithm makes sure that the number of considered boundary vectors in Eq. (6) becomes smaller with every iteration. The idea is, that the influence of distant vectors is negligible. Hence, they can be omitted completely.
Fig. 2 shows how the progressive refinement algorithm works. Only vectors from previous steps within the gray square region around a point P are included in the computation. Vectors to be computed in the current iteration are drawn in red and the exemplarily considered point P is marked by a red ring. The pre-computed boundary vectors are drawn in black. Vectors computed in earlier refinement steps are called raster vectors and are drawn in green. The set of raster vectors gradually becomes finer until the desired resolution is reached.
To get periodic boundary conditions, in the first iteration (Fig. 2a) the region of considered NVs covers the structure itself and its four neighboring structures. Only the boundary vectors are used as an input for the interpolation in this step, since no raster vectors exist, yet. In the second iteration (Fig. 2b) the region size is identical to the size of the unit cell. Now the red vectors are computed from the boundary vectors in the surrounding square and the raster vector (only a single one at this point). In the third iteration (Fig. 2c) the region shrinks to a quarter of the structure size, leaving only five boundary vectors to consider plus four raster vectors from the iteration before in the corners. Another quartering is done in the last iteration leaving only one boundary vector to consider plus the four raster vectors in the corners. (Fig. 2d). Of course, in practice the sampling of the initial boundary vectors and the final raster size are chosen to be finer, which means a much higher number of iteration steps. As mentioned this algorithm features a complexity of only 𝒪(R2) which gives the needed performance boost to calculate NV fields in a fraction of the time the RCWA takes.
5. Results
In the following we present the simulation results of two example structures. The NV-method presented here is compared to the zigzag method of Li [7] as well as to an implementation which ignores Li’s rules which we refer to as the “original” formulation. Our implementation is based on [13] which is a more modern formulation of the authors of [1] and at the same time one of the most appreciated publications on RCWA. Although 2D periodic gratings are not explicitly mentioned there, its consideration in RCWA is straightforward as pointed out first by Bräuer and Bryngdahl [14].
For both examples the simulation parameters from Table 1 are taken. The first example is a
Table 1. Simulation parameters
C-shaped structure. Fig. 3 shows the two variants with and without preferred direction of the NVs. Plotting the corresponding diffraction efficiencies as a function of the retained Fourier modes, commonly referred to as “truncation order” (Fig. 4), one can see, that the NV method has a considerably faster convergence than both the original formulation and the zigzag method. It also shows that both variants, although having a completely different look, yield almost the same convergence rate.
The second example is a structure which is built out of two intersecting circles. In contrast to the C-shaped structure which is “hand drafted” it features a certain regularity. However, in order to make sure that the method succeeds not only due to symmetry, odd values for the circle parameters were chosen. For convenience these data are reproduced in Table 2 with the origin being in the upper left corner of the unit cell. NV fields for this structure are shown in Fig. 5. Again, the NV method shows better convergence than the other methods (Fig. 6). Furthermore it is obvious that none of the two NV variants are superior to the other.
A general rule of thumb cannot be made for which of both variants must be used. It is also impossible to make rules by arguing with the number of singularities in the NV field or with the distance of singularities from boundaries. Instead, there is very strong evidence that the proper
Table 2. Data of the intersecting circle structure
representation of N2x, N2y and NxNy as truncated Fourier series plays an important role. To be precise, it is important that the deviations of the truncated Fourier series of N2x, N2y and NxNy from the exact values are as small as possible at the material boundary. A detailed investigation of this matter is presented in [15] but omitted here for the sake of brevity.
Fig. 3. NV field for the C-shaped structure (a) without preferred direction, (b) with preferred direction to the right. (Note that both vector fields have in fact a higher resolution than shown here)
Fig. 5. NV field for intersecting circles (a) without preferred direction, (b) with preferred direction
6. Conclusion
With the presented algorithm it is now possible to calculate NV fields for arbitrary structures, e.g. concave shapes or structures with several arbitrarily connected materials. While in this work we presented two variants of the algorithm, both yielding similar results, in practice one of these variants can be chosen. In conjunction with this algorithm the NV method with its excellent convergence is a very attractive alternative to Li’s zigzag method for crossed gratings.
References and links
1. M. G. Moharam and T. K. Gaylord, “Rigorous coupled-wave analysis of planar-grating diffraction,” J. Opt. Soc. Am. 71, 811–818 (1981). [CrossRef]
2. L. Li and C. W. Haggans, “Convergence of the coupled-wave method for metallic lamellar diffraction gratings,” J. Opt. Soc. Am. A 10, 1184–1189 (1993). [CrossRef]
3. P. Lalanne and M. Morris, “Highly improved convergence of the coupled-wave method for TM polarization,” J. Opt. Soc. Am. A 13, 779–784 (1996). [CrossRef]
4. G. Granet and B. Guizal, “Efficient implemenation of the coupled-wave method for metallic lamellar grating in TM polarization,” J. Opt. Soc. Am. A 13, 1019–1023 (1996). [CrossRef]
5. L. Li, “Use of Fourier series in the analysis of discontinuous periodic structures,” J. Opt. Soc. Am. A 13, 1870–1876 (1996). [CrossRef]
6. P. Lalanne, “Improved formulation of the coupled-wave method for two-dimensional gratings,” J. Opt. Soc. Am. A 14, 1592–1598 (1997). [CrossRef]
7. L. Li, “New formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. Soc. Am. A 14, 2758–2767 (1997). [CrossRef]
8. T. Schuster, J. Ruoff, N. Kerwien, S. Rafler, and W. Osten, “Normal vector method for convergence improvement using the RCWA for crossed gratings,” J. Opt. Soc. Am. A 24, 2880–2890 (2007). [CrossRef]
9. M. Nevière and E. Popov, Light Propagation in Periodic Media - Differential Theory and Design (Marcel Dekker, Inc., New York, 2003).
10. E. Popov and M. Nevière, “Maxwell equations in Fourier space: fast-converging formulation for diffraction by arbitrary shaped, periodic, anisotropic media,” J. Opt. Soc. Am. A 18, 2886–2894 (2001). [CrossRef]
11. D. Shepard, “A two-dimensional interpolation function for irregularly-spaced data,” pp. 517–524 (ACM National Conference, 1968).
12. R. Franke, “Scattered Data Interpolation: Tests of Some Methods,” Math. Comput. 38, 181–200 (1982).
13. M. G. Moharam, T. K. Gaylord, E. B. Grann, and D. A. Pommet, “Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings,” J. Opt. Soc. Am. A 12, 1068–1076 (1995). [CrossRef]
14. R. Bräuer and O. Bryngdahl, “Electromagnetic diffraction analysis of two-dimensional gratings” Opt. Commun. 100, 1–5 (1993). [CrossRef]
15. P. Götz, “Simulation von Lichtbeugung mit der Normalenvektormethode,” Diploma thesis, Universität Stuttgart, Germany (2008).
