Reformulation of the Fourier modal method with adaptive spatial resolution: application to multilevel profiles

T. Vallius; M. Honkanen

doi:10.1364/OE.10.000024

1 Introduction

Of the several rigorous methods frequently used to analyze periodic structures scattering light, the Fourier modal method (FMM) is perhaps the most versatile method available enabling reliable analysis of multilevel as well as continuous structures. Since the method was introduced by Knop [1], it has been applied to a variety of problems, in which rigorous analysis of micro structures has been necessary. A characteristic feature of the method is that it converges when the number of modes and Rayleigh orders in the analysis is increased yet requiring sometimes more modes than even up-to-date computers can efficiently calculate, being a considerable problem particularly in TM polarization.

The correct factorization rules of finite Laurent series [2, 3] derived by Li [4] improved the convergence in TM polarization, but still the problems due to the weak convergence of Fourier series with discontinuous permittivity distributions remained. The major disturbing factor was the Gibbs’ phenomenon arising at the discontinuous points damaging the convergence and it had to be diminished, or even avoided.

A parametric representation of the coordinate axis allowed a spatially adaptive resolution increasing the sampling in the neighborhood of the transitions and substantial improvements in the convergence were reported by Granet [5]. The method required solving eigenvalue problems in the homogeneous regions, i.e, in the regions before and after the modulated area, which enabled the matching of the boundary conditions in the adaptive coordinate system. The use of parametric representation also required that the resolution function had to be similar in every region; thus the transitions in multilevel cases had to be equal in every grating layer. However, with multilevel structures the spacing of the transitions in different grating layers may usually vary independently and similar representation of the parametric coordinates in every layer is not feasible. Therefore a different approach must be developed to enable the analysis of multilevel structures with spatially adaptive resolution.

We will modify the method introduced by Granet [5] and make a generalization to multilevel profiles with different transitions. Also the parametric representation is slightly modified for achieving better convergence of the eigenvalues. For the boundary problem modes projected on the Fourier basis in the x space will be used enabling different parametric representations in every modulated layer and the extra eigenvalue problems will be evaded. Projections of the field onto similar basis in every grating layer allows also arbitrary spacing of the transitions in each layer.

2 Theory

Let us assume the y-invariant geometry illustrated in fig. 1. The media in the regions 1 and 3 are assumed to be homogeneous with refractive indices n ₁ and n ₃, respectively. In region 2 the relative permittivity ∊(x) = n(x)² is assumed periodic in x-direction, with period d. Region 2 is also divided into J layers: in each layer the relative permittivity can be arbitrarily modulated in x-direction but is assumed z-invariant. The structure is illuminated by an infinite plane wave from the negative z-direction propagating in the xz-plane with an angle θ between the z-axis and the wave vector k, where |k| = k = 2π/λ, and λ is the wavelength in vacuum. Since the problem is completely y-invariant, all the partial y-derivatives vanish in Maxwell’s equations and one has the so-called TE/TM-polarization decomposition [6].

2.1 Field in unmodulated regions

Let us consider a field in the homogeneous regions 1 and 3. Since the refractive index distribution in region 2 is periodic the reflected and transmitted fields must be pseudo-periodic [6] and the fields can be expressed in the form

E_{1} (x, z) = \exp [i (α_{0} x + r_{0} z)] + \sum_{m} R_{m} \exp [i (α_{m} x - r_{m} z)],

E_{3} (x, z) = \sum_{m} T_{m} \exp [i (α_{m} x + t_{m} z)],

Fig. 1. y-invariant grating geometry.

Download Full Size | PDF

where

r_{m} = \sqrt{{(k n_{1})}^{2} - α_{m}^{2}},

t_{m} = \sqrt{{(k n_{3})}^{2} - α_{m}^{2}}

with α_m = α ₀ + 2πm/d, and α ₀ = kn ₁ sin θ. Corresponding representations for the magnetic field in TM polarization can be achieved by replacing E by H. R_m and T_m denote the complex amplitudes of the m’th reflected or transmitted diffraction order and the diffraction efficiencies are given by

η R_{m} = ℜ (r_{m} / r_{0}) {∣ R_{m} ∣}^{2}

η T_{m} = C ℜ (t_{m} / r_{0}) {∣ T_{m} ∣}^{2},

where ℜ denotes the real part. In TE polarization the parameter C = 1 and in TM polarization C = (n ₁/n ₃)².

2.2 Eigenvalue problem

Let us proceed by considering the field inside the modulated region 2. Wave equations for y invariant structures [1, 7] when light arrives in xz plane are

\frac{\partial}{\partial x} \frac{\partial E_{y}}{\partial x} + \frac{\partial}{\partial z} \frac{\partial E_{y}}{\partial z} + k^{2} n^{2} (x) E_{y} = 0

in TE polarization, and

n^{2} (x) {\frac{\partial}{\partial x} [\frac{1}{n^{2} x} \frac{\partial H_{y}}{\partial x}] + k^{2} H_{y}} + \frac{\partial}{\partial z} \frac{\partial H_{y}}{\partial z} = o

in TM polarization. These equations give us the eigenvalue problem used in the original FMM [1]. Since the permittivity distribution is discontinuous in the x-direction the Gibbs phenomenon appears and slow convergence of the Fourier series is observed. Thus better convergence rate can be achieved by choosing a new Fourier space, where the variable x is replaced by u. The dependence between the variables x and u is presented by parametric equations that lead to new eigenvalue equations.

We apply a change of variable x → x(u):

\frac{\partial}{\partial x} \to \frac{\partial u}{\partial x} \frac{\partial}{\partial u},

and denote the resolution function by

f (u) = \frac{\partial x}{\partial u} .

Two auxiliary functions are defined in the form [5]

a (u) = n^{2} (u) f (u),

b (u) = \frac{f (u)}{n^{2} (u)} .

Also in the u space the fields E_y and H_y are pseudoperiodic with period d. Hence each mode can presented in the form

E_{y} (u, z) = \exp (i γ z) \sum_{m} E_{m} \exp (i α_{m} u),

where E_m is an eigenvector containing the Fourier coefficients of the mode, and γ is an eigenvalue that defines the propagation in z direction. Similar representation is valid also for the modes of the magnetic field H_y.

We substitute the representation (13) to the wave equations (7)–(8), calculate the Fourier coefficients for the field applying the correct rules for the products of the series of discontinuous functions [4], and use the equations (9), (10), (11), and (12), whereupon we achieve the eigenvalue problems [5]

γ^{2} E = f^{- 1} [k^{2} a - α f^{- 1} α],

for TE polarization, and

γ^{2} H = b^{- 1} [k^{2} f - α a^{- 1} α],

for TM polarization, where f, a, and b are toeplitz matrices formed from the Fourier coefficients of the respective functions, and α is a diagonal matrix formed from α_m. The sign rule for the eigenvalues is known as ℜ{γ} + ℑ{γ} > 0, where ℑ denotes the imaginary part. The exact eigenvalues are independent of the chosen representation of the coordinates; the eigenvectors E and H depend on the coordinate system and must be transformed to a more convenient one.

2.3 Parametric representation of the coordinate x

We present the coordinate x as a function of u and the transition points are denoted by x_l in the x space and by u_l in the u space. Between the transitions l and l - 1 we use the function x_l(u) for the mapping between different spaces:

x_{l} (u) = a_{1} + a_{2} u + \frac{a_{3}}{2 π} \sin [2 π \frac{u - u_{l - 1}}{u_{l} - u_{l - 1}}],

where

a_{1} = \frac{u_{l} x_{l - 1} - u_{l - 1} x_{l}}{u_{l} - u_{l - 1}}

a_{2} = \frac{x_{l} - x_{l - 1}}{u_{l} - u_{l - 1}}

a_{3} = G (u_{l} - u_{l - 1}) - (x_{l} - x_{l - 1})

and G = f(u _l-1) = f(u_l). The difference between the representation (16) and the original parametric function [5, 8] is the modification of the term a ₃ allowing different spacing of the transitions in x and u spaces without discontinuities of the resolution function, which was our motivation for the reformulation of the method. In the original formulation the parameter a ₃ was defined as a ₃ = (G - 1)(x_l - x _l-1).

2.4 Boundary value problem

After solving the eigenvalue problems (14)–(15) in each layer the field in layer j can be written as a superposition of the modes given by Eq. 13:

E_{y} (u, z) = \sum_{v = 1}^{V} {A_{v, j} \exp [i γ_{v, j} (z - h_{j})] + B_{v, j} \exp [- i γ_{v, j} (z - h_{j + 1})]}

in TE polarization, and corresponding representation in TM polarization is achieved by replacing E by H. h_j’s are the height transitions in the z-direction. The field representations in each layer are in different bases exp[iα_m u_j(x)] which depend on the locations of the transitions in the respective layer. Thereby we expand the eigenfunction in each layer in terms of similar base functions in the x space. The boundary conditions between the layers j and j + 1 in TE polarization are

E_{y, j} (x, h_{j}) = E_{y, j + 1} (x, h_{j})

\frac{\partial}{\partial z} [E_{y, j} (x, h_{j})] = \frac{\partial}{\partial z} [E_{y, j + 1} (x, h_{j})]

and in TM polarization

H_{y, j} (x, h_{j}) = H_{y, j + 1} (x, h_{j})

\frac{\partial}{\partial z} [\frac{1}{∊_{j} (x)} H_{y, j} (x, h_{j})] = \frac{\partial}{\partial z} [\frac{1}{∊_{j + 1} (x)} H_{y, j + 1} (x, h_{j})] .

The orthogonal Fourier base in the x space is an apparent choice for the basis function allowing the use of general algorithms for calculating the projections, which are easily obtained by calculating the integral

K = {[K]}_{p, m} \frac{1}{d} \int_{0}^{d} f (u) \exp [- i α_{p} x (u) + i α_{m} u] d u,

and by multiplying the eigenvectors with the matrix K:

E_{j} = K E^{u}, H_{j} = K H^{u}, Q_{j} = K Q^{u} .

The matrices containing the eigenvectors in the x space are denoted by E _j and H _j in layer the j and the matrix Q _j includes the Fourier coefficients of the function Q_y,j(x) = H_y,j(x)/∊_j(x) in the equation (24); the superscript u refers to the eigenmatrices in the u space. By substituting the matrices in the x space to the boundary conditions (21)–(24) one obtains the following matrix equation in TM mode

[\begin{matrix} H_{j + 1} & - H_{j} \\ Q_{j + 1} Γ_{j + 1} & Q_{j} Γ_{j} \end{matrix}] [\begin{matrix} A_{j + 1} \\ B_{j} \end{matrix}] = [\begin{matrix} H_{j + 1} X_{j} & - H_{j} X_{j + 1} \\ Q_{j + 1} Γ_{j + 1} X_{j} & Q_{j} Γ_{j} X_{j + 1} \end{matrix}] [\begin{matrix} A_{j} \\ B_{j + 1} \end{matrix}],

Fig. 2. Convergence of the real part of the first eigenvalue in TE (a) and TM (b) polarization calculated by using the old FMM (---), the first parametric representation (⋯), and the new formulation (-).

Download Full Size | PDF

where Γ _j and X _j are diagonal matrices with elements γ_v,j and exp[iγ_v,j(h _j+i - h_j)], while the vectors A _j and B _j are defined in equation (20). The matrix equations for TE modes are obtained by replacing H and Q with E. For a numerically stable treatment of the evanescent waves at the boundaries we use the S matrix algorithm [9], which is probably the most beautiful and elegant way of solving the boundary value problem.

We stress that for the stability of the method the product Q(x) = H(x)/∊(x) must be calculated in the u space and transformed to the x space by Eq. (26), not by using the Laurent’s rule in the x space.

3 Numerical examples

In this section we investigate the convergence properties of the method compared to the original formulations of the FMM for multilevel structures [11, 7]. The convergence of eigenvalues as well as diffraction efficiencies will be studied and compared using different example structures.

3.1 Convergence of eigenvalues

The accuracy of the eigenvalues has a significant effect on the waveguide modes used in the boundary value problem defining the coupling efficiency between the modes in different grating layers. However, the convergence rate of the eigenvalues of gratings analyzed with the FMM is sometimes too slow when the accuracy of several decimals is needed.

Next we will proceed by considering a grating with period d = 1α containing a single transition at c_x = 0.9d. The parameter G is set to G = 0.001 for all the examples presented here. The refractive indices around the transition are n _2a = 1 and n _2b = 5+15i and normal incidence is assumed. We calculate the eigenvalues by using a growing number of the truncation order and compare the convergence of the real part of the first eigenvalue. In fig. 2a and 2b we see how the parametric representation of the x axis leads to faster convergence of the eigenvalue. We have also compared the convergence of the eigenvalue with two different locations of the transition in the u space. The dotted line illustrates the error when similar transitions in both x and u space have been used, whereas the solid line represents the convergence when the transition in the u space is located at c_u = 0.5d The use of equal transitions in each space gives us the parametric representation used in the original formulation [5]. Although the performance of the original parametric representation was good in ref. [5], it certainly has problems when the feature size is small, since they can not be scaled in the u space with differentiable mapping.

Fig. 3. Logarithmic decay of the error |γ ₃ - γ _exact| of the third eigenvalue in TE (a) and TM (b) polarization calculated by using the old FMM (- - -), the first parametric representation (⋯), and the new formulation (-).

Download Full Size | PDF

Fig. 4. The dependence between the coordinates x and u by using transition c_x = 0.9d in the x space and in the u space c_u = 0.5d which represents the new formulation (-) and c_u = 0.9 which represents the old formulation (---).

Download Full Size | PDF

Figure 3 shows the logarithmic error of the third eigenvalue of the grating as a function of the truncation order calculated by the Fourier modal method (dashed line) and by the new parametric formulation (solid line). Also the old parametric formulation [5] is used (dotted line) but it does not allow the spreading of the small features in the u space, like the reformulated representation does. Thus the Gibb’s phenomenon will ruin the shape of the small detail. We can observe that the new parametric representation converges faster to the reliable value without substantial fluctuations. The difference between the representations is illustrated figure 4 where both the new and the old formulations are presented. We stress that the reformulation of the parametric representation is not necessarily the optimal one since several other mappings than the aforementioned sine function might prove more suitable. This primarily serves as a new basis for research aimed at obtaining better convergence of the FMM.

The faster convergence rate takes place in the set of the eigenvalues that includes approximately third part of the calculated eigenvalues, as was noted also in the method introduced by Morf using orthogonal polynomials for solving diffraction problems [10]. Hence we neglect the eigenvalues whose convergence is not guaranteed using only one third of the values with smallest absolute values [10]. Thus, when we mention the truncation order, we mean the eigenvalue problem of the size 3M × 3M; the number of the Fourier coefficients and the modes is M, respectively.

Fig. 5. A multilevel conducting grating with parameters d = 1, h_j = 0.125_j, n_a = 0.1217 + 3.2966i, and n_b, = 1. Both n ₁ and n ₃ equal to 1.

Download Full Size | PDF

Fig. 6. The diffraction efficiency of the zeroth transmitted order in TE polarization with normal incidence for the grating in fig. 5 as a function of the wavelength α according to the FMM (a) and the parametric representation (b) with different number of the modes: 5 modes (-) and 10 modes (◇ ◇ ◇). The dashed line is the accurate result calculated by the FMM with 240 modes.

Download Full Size | PDF

Fig. 7. Same as fig. 6 in TM polarization. The number of the modes: 12 modes (-) and 24 modes (◇ ◇ ◇). The dashed line: the FMM with 240 modes.

Download Full Size | PDF

3.2 Accuracy of the diffraction efficiency

We consider a multilevel profile divided into 4 z invariant layers as illustrated in fig. 5 to determine convergence properties of diffraction orders. Results are presented in fig. 6 for TE polarized light, where the diffraction efficiency of the zeroth order analysed with FMM in fig. 6a is smaller than the correct value when five modes have been included in the analysis, but the method with parametric representation in fig. 6b is reliable even with 5 modes. The accurate results calculated by the FMM with 240 modes (---) agree well with the values of the parametric formulation. Figure 7 illustrates the corresponding results in TM polarization and remarkable changes in the curves with different number of modes can be observed in fig. 7a with FMM, while the parametric formulation lead to accurate values with smaller number of modes in fig. 7b.

Fig. 8. A checkerboard grating with parameters n_a = 5, n_b, = 1.5, h ₁ = 10, h ₂ = 20, and d = 1. The refractive indices in the homogeneous regions are n ₁ = 1 and n ₃ = 1.

Download Full Size | PDF

Fig. 9. The diffraction efficiency of the zeroth transmitted diffraction order of the checkerboard grating in TM polarization as a function of the wavelength with normal incidence. The curves have been calculated by the FMM (a) and the parametric representation (b) with19 (-) and 38 ( opex-10-1-24-i001 ) modes.

Download Full Size | PDF

In the second example we analyze a checkerboard grating profile shown in the fig. 8. This kind of structures cause sharp resonance peaks for certain wavelengths the accurate locations of which may be difficult to define. Figure 9 illustrates the efficiency of the zeroth diffraction order as a function of the wavelength using two different number of the modes in TM polarization. As was noticed also in ref. [10], the traditional multilevel FMM is not very reliable with this kind of profile because even small errors in the eigenvalues cause shifts of the resonance wavelengths when the number of the modes is increased (figure 9a). Nevertheless, the formulation using parametric representation of the coordinates is stable illustrating the spectral resolution reliably in fig. 9b. Especially the resonances at the wavelengths λ = 1.171 and λ = 1.161 have been shifted in the analysis with the FMM in fig. 9a, whereas in fig. 9b the locations remain the same.

In our last example case we consider a grating formed from cylinders in a high refractive material (Figure 11) for TE polarized plane wave under normal incidence. The grating structure has been divided into 120 slices of equal size and the refractive indices of the cylinders and the surrounding medium are n_a = 1 and n_b = 5. Figure 11 illustrates the efficiency of the zeroth transmitted diffraction order analyzed with FMM in fig. 11a and with the present method in fig. 11b. Note the evident change of the resonance peaks with FMM which is not observable when the adaptive resolution has been applied.

Fig. 10. A cylindrical grating with parameters n ₁ = n ₃ = 1, n_a = 1, n_b, = 5, h ₁ = 0.5, h ₂ = 1,d=1, and r = 0.25.

Download Full Size | PDF

Fig. 11. The diffraction efficiency of the zeroth transmitted diffraction order of the cylindrical grating in TE polarization as a function of the wavelength with normal incidence. The curves have been calculated by the FMM (a) and the parametric representation (b) with7 (-) and 14 (() modes.

Download Full Size | PDF

4 Conclusions

We have utilized the parametric representation of the coordinate system to enhance the reliability of the Fourier modal method for multilevel gratings. The mapping between the x and u spaces was slightly reformulated to enable different spacing of the transitions in each space. The eigenvalue problem was solved in the adaptive u space and the eigenvectors were transformed into the x space for solving the boundary conditions between different grating layers. Good convergence rates both in TE and in TM polarizations were achieved and the method proved reliable also for structures causing strong resonance peaks.

The method resembles the Morf’s method [10] because the size of the eigenvalue problem is unfortunately three times larger and the convergence rate of the eigenvalues is of the same kind. The calculation time is longer, but the accurate eigenvalues are close to the exact ones giving better reliability with structures that are difficult to analyze with the traditional FMM due to sharp resonance peaks. Thus the method would be useful for analyzing structures like photonic crystals that characteristically produce resonances with narrow wavelength bands.

References and links

1. K. Knop, “Rigorous diffraction theory for transmission phase gratings with deep rectangular grooves,” J. Opt. Soc. Am. 68, 1206–1210 (1978). [CrossRef]

2. P. Lalanne and G. M. Morris, “Highly improved convergence of the coupled-wave method for TM polarization,” J. Opt. Soc. Am. A 13, 779–784 (1996). [CrossRef]

3. G. Granet and B. Guizal, “Efficient implementation for the coupled-wave method for metallic lamellar gratings in TM polarization,” J. Opt. Soc. Am. A 13, 1019–1023 (1996). [CrossRef]

4. L. Li, “Use of Fourier series in the analysis of discontinuous periodic structures,” J. Opt. Soc. Am. A 13, 1870–1876 (1996). [CrossRef]

5. G. Granet, “Reformulation of the lamellar grating problem through the concept of adaptive spatial resolution,” J. Opt. Soc. Am. A 16, 2510–2516 (1999). [CrossRef]

6. R. Petit, ed., Electromagnetic theory of gratings (Springer-Verlag, Berlin, 1980). [CrossRef]

7. J. Turunen, “Diffraction theory of microrelief gratings,” Chap 2 in Micro-Optics: Elements, Systems and Applications,H. P. Herzig, ed. (Taylor & Francis, Cornwall, 1997)

08. G. Granet, J. Chandezon, J.-P. Plumey, and K. Raniriharinosy, “Reformulation of the coordinate transformation method through the concept of adaptive spatial resolution. Application to trapezoidal gratings,” J. Opt. Soc. Am. A 18, 2102–2108 (2001). [CrossRef]

9. L. Li, “Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings,” J. Opt. Soc. Am. A 13, 1024–1035 (1996). [CrossRef]

10. R. H. Morf, “Exponentially convergent and numerically efficient solution of Maxwell’s equations for lamellar gratings,” J. Opt. Soc. Am. A 12, 1043–1056 (1995). [CrossRef]

11. D. Nyyssonen and C. P. Kirk, “Optical microscope imaging of lines patterned in thick layers with variable edge geometry,” J. Opt. Soc. Am. A 5, 1270–1280 (1988). [CrossRef]

Reformulation of the Fourier modal method with adaptive spatial resolution: application to multilevel profiles

Abstract

1 Introduction

2 Theory

2.1 Field in unmodulated regions

2.2 Eigenvalue problem

2.3 Parametric representation of the coordinate x

2.4 Boundary value problem

3 Numerical examples

3.1 Convergence of eigenvalues

3.2 Accuracy of the diffraction efficiency

4 Conclusions

References and links

Cited By

Figures (11)

Equations (28)

Optics Express