## A construction guide to analytically generated meshes for the Fourier Modal Method |

Optics Express, Vol. 20, Issue 16, pp. 17319-17347 (2012)

http://dx.doi.org/10.1364/OE.20.017319

Acrobat PDF (12084 KB)

### Abstract

The concepts of adaptive coordinates and adaptive spatial resolution significantly enhance the performance of Fourier Modal Method for the simulation of periodic photonic structures, especially metallo-dielectric systems. We present several approaches for constructing different types of analytical coordinate transformations that are applicable to a great variety of structures. In addition, we analyze these meshes with an emphasis on the resulting convergence characteristics. This allows us to formulate general guidelines for the choice of mesh type and mesh parameters.

© 2012 OSA

## 1. Introduction

1. K. Busch, G. von Freymann, S. Linden, S. F. Mingaleev, L. Tkeshelashvili, and M. Wegener, “Periodic nanostructures for photonics,” Phys. Rep. **444**, 101–202 (2007). [CrossRef]

*z*-direction (propagation direction) and infinitely extended into the

*xy*-plane, where the dielectric material properties are lattice-periodic with respect to the lateral coordinates in the

*xy*-plane. The basic idea of FMM is that this periodicity should be reflected in the basis functions that are used to solve Maxwell’s equations - an appropriately constructed plane wave basis is the natural choice for representing the discrete periodicity with respect to the lateral coordinates. The three-dimensional system is stair-cased in propagation direction, i.e., it is sliced into several two-dimensional layers each of which is assumed to be invariant in the propagation direction. This invariance allows for a harmonic ansatz for the

*z*-dependence of the electromagnetic field, thus facilitating a reformulation of Maxwell’s equations into an eigenvalue problem for the corresponding propagation constant. After solving the corresponding eigenvalue problems, the electromagnetic fields in the different layers are expanded into the respective eigenmodes and the layers are recursively connected with a scattering matrix algorithm, which ensures the continuity conditions of the tangential fields in adjacent layers [2

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

3. 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]

4. G. Granet and B. Guizal, “Efficient implementation of the coupled-wave method for metallic lamellar gratings 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]

10. S. Essig and K. Busch, “Generation of adaptive coordinates and their use in the Fourier Modal Method,” Opt. Express **18**, 23258–23274 (2010). [CrossRef] [PubMed]

6. T. Weiss, G. Granet, N. A. Gippius, S. G. Tikhodeev, and H. Giessen, “Matched coordinates and adaptive spatial resolution in the Fourier modal method,” Opt. Express **17**, 8051–8061 (2009). [CrossRef] [PubMed]

## 2. Covariant formulation of the Fourier Modal Method with generalized coordinates

6. T. Weiss, G. Granet, N. A. Gippius, S. G. Tikhodeev, and H. Giessen, “Matched coordinates and adaptive spatial resolution in the Fourier modal method,” Opt. Express **17**, 8051–8061 (2009). [CrossRef] [PubMed]

10. S. Essig and K. Busch, “Generation of adaptive coordinates and their use in the Fourier Modal Method,” Opt. Express **18**, 23258–23274 (2010). [CrossRef] [PubMed]

11. L. Li, “Fourier modal method for crossed anisotropic gratings with arbitrary permittivity and permeability tensors,” J. Opt. A **5**, 345–355 (2003). [CrossRef]

*Ox̄*

^{1}

*x̄*

^{2}

*x̄*

^{3}and a curvilinear coordinate system

*Ox*

^{1}

*x*

^{2}

*x*

^{3}. They are connected by a transformation within the lateral plane of the general form Our aim is to solve Maxwell’s curl equations which are given in covariant form as [11

11. L. Li, “Fourier modal method for crossed anisotropic gratings with arbitrary permittivity and permeability tensors,” J. Opt. A **5**, 345–355 (2003). [CrossRef]

*ξ*denotes the Levi-Civita symbol and

*E*and

_{σ}*H*are covariant components of the electric and magnetic field. Here, and throughout the manuscript, Greek indices run over the range 1, 2, 3, and we utilize the Einstein sum convention where repeated indices are implicitly summed over. Further

_{σ}*k*

_{0}=

*ω*/

*c*is the vacuum wave number (

*c*denotes the vacuum speed of light and

*ω*the frequency). The metric tensor is and

*g*(as used in Eqs. (4) and (5)) denotes the reciprocal of its determinant. As a consequence of the coordinate transformation, the permittivity changes according to where

*ε̄*is the permittivity tensor in the Cartesian system. The permeability transforms identically. At this point, we wish to emphasize that an isotropic material will in general become anisotropic after applying a coordinate transformation.

^{τχ}*x̄*

^{3}direction, and formulate an eigenvalue problem using the correct Fourier factorization rules [5

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

11. L. Li, “Fourier modal method for crossed anisotropic gratings with arbitrary permittivity and permeability tensors,” J. Opt. A **5**, 345–355 (2003). [CrossRef]

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

12. L. Li, “New formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. Soc. Am. A **14**, 2758–2767 (1997). [CrossRef]

10. S. Essig and K. Busch, “Generation of adaptive coordinates and their use in the Fourier Modal Method,” Opt. Express **18**, 23258–23274 (2010). [CrossRef] [PubMed]

**5**, 345–355 (2003). [CrossRef]

13. 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]

*ε*

_{fib}= 2 in a background material with

*ε*= 1 on a grid with 1024 × 1024 points. This is the same discretization we use in our computations in section 9. A coordinate transformation designed for a cylinder is given in Ref. [6

_{bg}6. T. Weiss, G. Granet, N. A. Gippius, S. G. Tikhodeev, and H. Giessen, “Matched coordinates and adaptive spatial resolution in the Fourier modal method,” Opt. Express **17**, 8051–8061 (2009). [CrossRef] [PubMed]

*x̄*

^{1}- and in

*x̄*

^{2}-direction are parallel. In the following, we show how this transformation affects the eigenvalue problem of the FMM. The matrix that is Fourier transformed when using FMM with AC and/or ASR is not the permittivity itself but rather

*ε*like in Eq. (7). Therefore, we refer to

^{ρσ}*effective permittivity*. For our present example, we depict the

*∂x*

^{1}/

*∂x̄*

^{1}. It can be obtained by (for other derivatives, the situation is completely analogous) If the coordinate lines of the mesh are parallel in

*x̄*

^{1}- and

*x̄*

^{2}-direction at one point, the determinant of the Jacobian

*J*of the coordinate transformation will become zero at this point. This leads to a diverging derivative

*∂x*

^{1}/

*∂x̄*

^{1}and, thereby, to a diverging effective permittivity.

*nondifferentiable meshes*even for complex structures in section 3. The second class of meshes, called

*smoothed meshes*, is built in such a way that the four specific coordinate lines are smooth. In section 4, we discuss that this does not yet lead to fully differentiable meshes. However, smoothed meshes deform the shape of the structure in the transformed space away from the advantageous rectangular forms. Finally, we construct a third class of meshes for which all coordinate lines are differentiable everywhere in the unit cell. This

*differentiable meshing*is discussed in section 5.

## 3. Construction of nondifferentiable meshes

*r*

_{1}and the center (

*M*

_{1},

*N*

_{1}) of the large (outer) circle, the radius

*r*

_{2}and the center (

*M*

_{2},

*N*

_{2}) of the small (inner) circle, and the radii

*r*

_{3}and

*r*

_{4}of the tiny circles that form the tips (see Fig. 3(b)). We assume that these quantities are given and, furthermore, that the crescent is symmetric (i.e.,

*r*

_{4}=

*r*

_{3}) and not rotated in the unit cell, i.e.,

*M*

_{1}=

*M*

_{2}. While this makes the following example easier, it does in no way restrict the method. Henceforth, all geometric quantities such as radii, side lengths, etc., are given in units of the lattice constant.

*x*∈ [

*c*,

*d*], defines a straight line through the points (

*c*,

*c̄*) and (

*d*,

*d̄*). Moreover, we can enter entire coordinate lines in

*LT*(i.e.,

*c*=

*c*(

*y*), etc.) such that

*LT*describes the linear transition between these coordinate lines. Therefore, when a coordinate line

*c*(

*y*) is mapped on a new coordinate line

*c̄*(

*y*) and a coordinate line

*d*(

*y*) is mapped on a new coordinate line

*d̄*(

*y*) then the mapping of the coordinate line

*x*is given by

*LT*. In Fig. 4, we provide an illustration of such a mapping where the left and right vertical coordinate lines are mapped onto curvy lines (green and blue). The coordinate lines in between (marked with crosses) are given by a linear transition from the left to the right curvy line, i.e., the distances on the horizontal axis are the same between the coordinate lines. As in any linear function, it does not matter in which sequence the points (or coordinate lines) are arranged—i.e.,

*LT*(

*c*,

*c̄*,

*d*,

*d̄*,

*e*) ≡

*LT*(

*d*,

*d̄*,

*c*,

*c̄*,

*e*).

*P*,

*Q*,

*R*and

*S*are found. Point

*P*is given by Eq. (11) and point

*Q*and

*R*follow from the structure’s symmetry. The (Cartesian) coordinate lines passing through any two of these characteristic points divide the unit cell into the zones ① to ⑨. Second, the so selected (Cartesian) coordinate lines have to be deformed in order to follow the crescent’s surface and this surface has to be parametrized (see Eq. (9)). Third and final, we have to find the two-dimensional transformation function (

*x̄*

^{1}(

*x*

^{1},

*x*

^{2}),

*x̄*

^{2}(

*x*

^{1},

*x*

^{2})) that maps the unit cell onto itself by a linear transition between these deformed coordinate lines. This is done separately for all nine zones in the unit cell and we will discuss only zone ④ (see Figs. 5(a) and 5(b)). The red (horizontal) coordinate lines in zone ④ in Fig. 5(a) are not affected by the mapping (i.e. they are in the same position in Figs. 5(a) and 5(b)). This means that the mapping for

*x̄*

^{2}(corresponding to the horizontal red lines) is the identity mapping The same applies to the leftmost coordinate line. The rightmost coordinate line experiences a deformation from a straight line to the arc given in Eq. (9). Explicitly, this reads

*LT*as described above. Analogous mappings can be constructed for all other zones in the unit cell [0, 1] × [0,

*L*].

*L*,

*R*,

*T*, and

*B*refer to left, right, top, and bottom, respectively.

*x̄*

^{2}mapping is more complex but the construction follows the same principles. It reads

**17**, 8051–8061 (2009). [CrossRef] [PubMed]

## 4. Construction of smoothed meshes

*x̄*

^{1}- and in

*x̄*

^{2}-direction (see Fig. 2(b)). The latter causes a divergent effective permittivity (see Eq. (7)). Therefore, we will now describe a way to design smoothed meshes that avoid these singularities. In this context,

*smoothed*means that the characteristic coordinate lines are mapped onto smooth curves. The result are meshes whose partial derivatives are still discontinuous at several points. In this sense, the present section represents an intermediate step from nondifferentiable to fully differentiable meshes which are discussed in the subsequent section.

*r*located in the center. Just as in the previous section, we choose characteristic points, namely the four points

*R*= (

*x*

_{−},

*x*

_{−}),

*S*= (

*x*,

_{+}*x*

_{−}),

*P*= (

*x*

_{−},

*x*

_{+}) and

*Q*= (

*x*

_{+},

*x*

_{+}), where we have introduced the abbreviations

*a*and at point

*ū*. We determine

*ū*by

*choosing*a smoothing parameter

*τ*, i.e.,

*ū*=

*x*

_{−}+

*τ*. Apparently, this adjustment of the transition function is not unique and one could use functions other than parabolas. By assuming a general form of for the parabola

*g*we assure continuity and differentiability at

*a*. Demanding continuity and differentiability at

*ū*determines the parameters

*a*and

*b*. The mapping itself is performed with the same methods as described in the preceding section. The coordinate lines that pass through any pair of characteristic points

*P*,

*Q*,

*R*,

*S*are mapped onto the smoothed curves as depicted in Fig. 8. The other coordinate lines follow by a linear transition. The functions needed are given in Eq. (13), Eq. (20) and by

*x̄*

^{2}then reads

*x̄*

^{1}can be found in an analogous manner. For our simple system, symmetry mandates that the coordinate

*x̄*

^{1}is constructed from the horizontal line mapping, Eq. (23), by interchanging

*x*

^{1}and

*x*

^{2}. We depict in Fig. 9 the final results of this mapping (cf. Fig. 2(b) for a corresponding ”unsmoothed” mesh). An important effect of this smoothing procedure is that the structure’s surface is no longer grid-aligned in the undistorted, Cartesian space described by Eq. (7).

*x*

^{1}∈ (

*a*, 1 −

*a*),

*x*

^{2}=

*x*

_{±}and

*x*

^{1}=

*x*

_{±},

*x*

^{2}∈ (

*a*, 1 −

*a*). The effect of this nondifferentiability manifests itself in a sudden change of the density of coordinate lines (see Fig. 9). As mentioned above, without any smoothing (i.e., Eq. (23) with

*τ*= 0) we obtain the expressions for the nondifferentiable circle mesh as described in Ref. [6

**17**, 8051–8061 (2009). [CrossRef] [PubMed]

*nondifferentiable circle mesh*. We will now improve the smoothed meshes by illustrating how to construct fully 2D differentiable meshes, i.e., meshes for which all partial derivatives exist and are continuous at all points throughout the unit cell.

## 5. Construction of differentiable meshes

*all*coordinate lines. We reuse the illustrative example of the previous section: A circle with radius

*r*is centered within a square unit cell [0, 1] × [0, 1]. The intersection of parabola and straight line is given by

*a*just like in Eq. (21). Also, we use the same definition of

*x*

_{±}and of

*ū*=

*x*

_{−}+

*τ*. Again,

*τ*is the only free parameter and determines how strongly the mesh is smoothed. In Fig. 10(a), we display how the unit cell is divided for the construction of the mapping. We restrict our discussion to the mapping in the zones ① to ⑥ —the other regions may be treated in an analogous fashion or follow directly from symmetry considerations.

**Zones ①, ②, ③:**In these regions, we choose the mappings exactly as in the case of the nondifferentiable circle mesh (see Ref. [6

**17**, 8051–8061 (2009). [CrossRef] [PubMed]

**Zone ④ :**To construct the mapping in this zone we have to deal with transitions from straight lines (zone ①) to ellipse arcs (zone ②). The general form of the ellipse arcs is given by with the parameters

*δ*and

*γ*. From this, we can infer the ellipse parameters in zone ② by comparing with the mapping in Eq. (25):

*γ*= 1. We aim to connect the straight lines (i.e., the identity mapping) in zone ① with the ellipse arcs by a family of parabolas of the form Note, that even though

*a*is still given by Eq. (21),

*b*is not a constant anymore as it has been in Eq. (22)—thus

*b*(

*x*

^{1}) parametrizes the family of parabolas. Continuity and differentiability at

*x̄*

^{2}=

*ū*require that

**Zone ⑤:**To find suitable differentiable mappings in the zones ⑤ and ⑥ we will utilize an alternative approach: We formulate the requirements in terms of continuity and differentiability and then make an ansatz to find appropriate functions. We start by considering the boundary conditions for the

*x̄*

^{1}(

*x*

^{1},

*x*

^{2}) component in zone ⑥ which follow directly from the transformations in the neighboring zones ② and ③. Continuity requires whereas differentiability requires

*f*which is assumed to be differentiable but otherwise still unknown. The ansatz reads This ansatz fulfills the continuity requirements of Eq. (30) provided

*f*does not exhibit a pole at

*x*

^{1}=

*a*or

*x*

^{1}=

*ū*. Further, the ansatz ensures that the derivative with respect to

*x*

^{2}fulfills Eq. (31) provided (

*∂f*/

*∂x*

^{2}) does not exhibit poles at

*x*

^{1}=

*a*or

*x*

^{1}=

*ū*. The remaining requirements refer to the derivative of Eq. (32) with respect to

*x*

^{1}which reads

*ū*and

*a*and comparing with Eq. (31) results in

*f*that fulfills Eq. (34) and does not have poles at

*x*

^{1}=

*a*or

*x*

^{1}=

*ū*is a suitable function so that the ansatz made in Eq. (32) fulfills all requirements in Eqs. (30) and (31). We choose a linear function which is given by Thereby, we acquired a differentiable mapping

*x̄*

^{1}(

*x*

^{1},

*x*

^{2}) in zone ⑤—namely Eq. (32) with

*f*given by Eq. (35).

*x̄*

^{2}(

*x*

^{1},

*x*

^{2}) in this zone by a ”parabola-construction” as in the case of zone ④, this time, however, in the

*x*

^{2}-direction Then,

*b*(

*x*

^{2}) is identical to Eq. (29) except that

*δ*is interchanged with

*γ*and

*x̄*

^{2}transformation in zone ③ (see Eq. (26)) are

**Zone ⑥:**In this zone we only need to find the mapping of one component because any point in this area obeys the symmetry relation The boundary conditions for this zone are derived from the mappings in zones ④ and ⑤. In particular, we use symmetry to connect the

*x̄*

^{2}mapping in zone ⑤ (Eq. (36)) to the

*x̄*

^{1}mapping in the area [

*ū*, 1−

*ū*]×[

*a*,

*ū*]. With the abbreviations Δ :=

*x*

_{+}−

*x*

_{−},

*v̄*:=

*ū*−0.5 and

*τ*. In Fig. 11 we depict the resulting meshes for two different values of

*τ*.

*x*

_{−}(see Fig. 8). In order to find the corresponding coordinate line

*α*for the differentiable mesh, we numerically compute the point that fulfills

## 6. Enforcing periodicity

*mirror structure*. We illustrate this approach in Fig. 14(c). Here, a mesh for the cross section of a trapezoidal rib waveguide has been created. The periodicity is assured by an additional structure in the unit cell which features similar geometrical properties as the structure under investigation and, thus, periodicity is restored. As this mirror structure has the same permittivity as the background medium it is physically nonexistent and the numerical artifacts created by the distorted mesh are negligible.

*μ*m × 4

*μ*m and the symmetric trapezoid has a top length of 240 nm and a bottom length of 960 nm. Its height is 360 nm. This trapezoid consists of silicon with a dielectric constant

*ε*= 12.1 while the background medium is air with

_{s}*ε*

_{bg}= 1. The incoming wavelength was chosen to be the telecom wavelength at 1550 nm.

## 7. Nonrectangular unit cells

12. L. Li, “New formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. Soc. Am. A **14**, 2758–2767 (1997). [CrossRef]

*α*and a rectangular unit cell is given by [12

12. L. Li, “New formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. Soc. Am. A **14**, 2758–2767 (1997). [CrossRef]

*x*,

_{α}*y*, and

_{α}*x*

_{rec},

*y*

_{rec}refer to the coordinates in ordinary (Cartesian) and transformed space, respectively. We illustrate this in Figs. 16(a) and 16(b), where, in addition, we also depict a circle which we would like to mesh. Upon mapping the nonrectangular unit cell onto a rectangular one, the circle that is centered in the unit cell turns into an ellipse. We can mesh this ellipse by means of any of the methods discussed in sections 3–5 (where applicable together with the methods to recover periodic meshes introduced in section 6). Mapping back onto ordinary (Cartesian) space then delivers the final mesh as displayed in Fig. 16(c). Actually, we have obtained this mesh by first compressing the coordinate lines (see section 8) of the ellipse mesh displayed in Fig. 14(b) and mapping back this mesh to ordinary (Cartesian) space.

*x*

_{0},

*y*

_{0}) that is rotated by an angle

*ϕ*with semi-axes

*c*and

*d*is described by Comparing Eq. (44) with Eq. (45) yields By symmetry, another solution which yields the same ellipse is

*ϕ*= 45° and

*c*and

*d*interchanged. The resulting ellipse is shown in Fig. 16(b). We have obtained the resulting mesh displayed in Fig. 16(c) using the mapping in Eq. (42) with

*α*= 30° (hexagonal lattice).

## 8. Compression of coordinate lines - adaptive spatial resolution

9. T. Vallius and M. Honkanen, “Reformulation of the Fourier modal method with adaptive spatial resolution: application to multilevel profiles,” Opt. Express **10**, 24–34 (2002). [PubMed]

9. T. Vallius and M. Honkanen, “Reformulation of the Fourier modal method with adaptive spatial resolution: application to multilevel profiles,” Opt. Express **10**, 24–34 (2002). [PubMed]

9. T. Vallius and M. Honkanen, “Reformulation of the Fourier modal method with adaptive spatial resolution: application to multilevel profiles,” Opt. Express **10**, 24–34 (2002). [PubMed]

*x*:

_{l}*l*= 1,...,

*n*} and {

*x̄*:

_{l}*l*= 1,...,

*n*} denote two sets of points. By design, this means that

*x̄*(

*x*) =

_{l}*x̄*for

_{l}*l*= 1,...,

*n*. Further,

*G*denotes the slope of the function at all positions

*x*. The compressed coordinates are described by

_{l}*x̄*(

*x*)—the transformation is a linear function superimposed with a sine. In the vicinity of the points

*x̄*the coordinate line density is increased—for instance, these points could represent material interfaces along the

_{l}*x*-coordinate line. The interval [

*x*

_{l}_{−1},

*x*] is mapped onto the interval [

_{l}*x̄*

_{l}_{−1},

*x̄*]. In order to maintain periodicity with period

_{l}*L*, one chooses

*x̄*

_{1}= 0 and

*x̄*=

_{n}*L*. This means that the coordinate lines are compressed at the unit cell edges. For the structures considered in Ref. [9

**10**, 24–34 (2002). [PubMed]

*x̄*

_{1}and

*x̄*

_{2}. Furthermore, we assume that the length of the interval that is mapped onto [

*x̄*

_{1},

*x̄*

_{2}] is given—it is referred to as Δ

*x*. Then, we perform the mapping Eq. (47) as if the first material interface lies at 0. Here, we have introduced the abbreviations The result is depicted in Fig. 17(a)—we now have a compression at two interfaces with the correct spacing between them, albeit at the wrong locations. In order to achieve the compression in the interior in the unit cell, we first shift the function by

*x̄*

_{1}upwards, see Fig. 17(b). The point where the compression function intersects with the unit cell edge, denoted

*x*̃, is Unfortunately, Eq. (52) is a transcendental equation and, therefore, we have to find the point

*x*̃ numerically. Next, we shift the function to the right by

*L*̃ =

*L*−

*x*̃, see Fig. 17(c). The part that is still outside the unit cell (marked in red) is attached at the left side of the unit cell—this makes sense due to periodic boundaries in the FMM. Finally, in Fig. 17(d) we depict how the equidistant, vertical coordinate lines are now mapped to the compressed horizontal lines. The compression function for two material interface points

*x̄*

_{1}and

*x̄*

_{2}in the interior of the unit cell then reads We want to stress that one needs to start with the separation Δ

*x*of the interfaces instead of starting with their exact locations

*x*

_{1}and

*x*

_{2}(as has been the case in Eq. (47)). This is due to the fact that these locations are determined by the shift of the entire function and cannot be specified beforehand.

*x̄*is monotonically increasing. If a constant coordinate line density is desired, the transformations above offer a simple way to do so. For two compression points,

*G*is chosen to be (

*x̄*

_{2}−

*x̄*

_{1})/Δ

*x*which leads to a vanishing prefactor of the sine in Eq. (47). This is how we have obtained the compression in Fig. 18(b).

## 9. Suitable parameter choice and convergence characteristics

*G*and Δ

*x*), we can use different types of meshes (Cartesian, nondifferentiable, smoothed, and differentiable) some of which have a smoothing parameter

*τ*, and we have to choose how many points are used for the real space discretization. As it is impossible to display every combination of those parameters, we aim at presenting the overall tendencies along with corresponding rules-of-thumb that we found in our extensive parameter scans.

### 9.1. Dielectric structures

*r*= 800 nm that is centered within a square unit cell and whose index will be varied. The lattice constant is

*d*= 4000 nm with a background material of

*ε*= 1. The number of plane waves will be varied, but they will always come from within a circle in reciprocal space centered around the origin. We analyze this system for a wavelength of

_{bg}*λ*= 800 nm.

*γ*is given by

*n*

_{eff}=

*λγ*/(2

*π*). First, we compare the different types of meshes for a small dielectric contrast. We initially refrain from compressing the coordinate lines in order to separate ASR and AC effects. In Fig. 19(a), we depict the convergence characteristics for a fiber with dielectric constant

*ε*

_{fib}= 2. The smoothed and especially the differentiable mesh performs better than the Cartesian or the nondifferentiable mesh. From now on, we will concentrate on the comparison of nondifferentiable and differentiable meshes since the smoothed meshes show a worse quality than the differentiable meshes when a compression is used.

*x̄*

^{1}- and

*x̄*

^{2}-direction. The material surface, i.e., the compression points are

*x̄*

_{1}=

*x*

_{−}and

*x̄*

_{2}=

*x*

_{+}for the nondifferentiable mesh and

*x̄*

^{1}=

*α*and

*x̄*

^{2}= 1 −

*α*for the differentiable mesh (see section 4 and 5, respectively). In Figs. 19(b) to 19(d) we display the dependence of the error on the compression parameters

*G*and Δ

*x*for the nondifferentiable mesh (

*τ*= 0, panel (b)) and the differentiable mesh with parameters

*τ*= 0.002 (panel (c)) and

*τ*= 0.015 (panel (d)). Each computation has been performed with 997 plane waves in conjunction with a discretization in transformed space of 1024 × 1024 points. In these graphs the parameter

*G*has been stepped in 0.01 intervals and Δ

*x*in 0.0025 intervals. We have used the same parameter stepping for

*G*and Δ

*x*in Figs. 19 to 22. For illustrative purposes, we display in Fig. 19(d) selected compression functions for vastly different parameter values

*G*and Δ

*x*.

*G*. This is intuitive since a small value of

*G*means that the coordinate line density at the surface of the structure is increased. However, we obtain the best results for this mesh for larger values of

*G*at specific values of Δ

*x*. For most values of these optimal combinations

*G*and Δ

*x*, a small change in Δ

*x*quickly compromises this optimal performance and significantly larger errors are obtained. This is in contrast to the characteristics of the differentiable mesh. For this mesh with

*τ*= 0.002, we infer from Fig. 19(c) very good performance for nearly all combinations of

*G*and Δ

*x*. The differentiable mesh with

*τ*= 0.015 exhibits a slightly worse performance but the qualitative behavior that very good results can be found for a wide parameter range is retained (see Fig. 19(d)). Another feature which we extract from Figs. 19(b) to 19(d) is that, for a fixed value of

*G*, the error exhibits an apparent oscillatory dependence on Δ

*x*. We will return to this issue once we have addressed the role of certain parameters.

*τ*= 0.002 (panel (b)) and with

*τ*= 0.015 (panel (c)). We do find our above (partial) summary confirmed: While the nondifferentiable mesh performs best it does so only for very specific combinations of

*G*and Δ

*x*, the differentiable meshes yield slightly worse results but for a much wider range of mesh parameters.

*ε*

_{fib}= 10. The background material is again

*ε*= 1. We have further checked that the above (partial) summary holds true for other values of the permittivity by performing calculations with varying

_{bg}*ε*

_{fib}where the compression parameters have been kept the same (not shown) as well as for varying wavelengths where all other parameters have been kept the same (not shown).

*τ*= 0.002 and

*τ*= 0.015 in Figs. 22(b) and 22(c), respectively, do not change much relative to Figs. 19(c) and 19(d). Of course, this has been expected as otherwise FMM computations would be rather useless to begin with. Nevertheless, throughout Fig. 22 we observe the same interesting feature which we have already pointed out in the context of Fig. 19: For fixed value of

*G*, the error apparently oscillates as a function of Δ

*x*. Therefore, we now turn to the analysis of this issue and investigate in more detail the parameter region delineated by the white box in Fig. 22(a), i.e., with a much higher resolution in

*G*and Δ

*x*.

*G*, it exhibits oscillations and even jumps as a function of Δ

*x*. This indicates that the error can be linked to the number of coordinate lines within the structure—recall, that Δ

*x*determines how many coordinate lines are inside the structure, while

*G*controls the density at the surface. In order to investigate this, we show in Fig. 23(b) the dependence on Δ

*x*of the relative error of the effective refractive index for each of the first ten eigenmodes for a fixed value

*G*= 0.165. In addition, we depict the distance between the physical surface and the numerical surface. The latter is given by the center of the two coordinate lines closest to the physical surface (for an illustration, see Fig. 23(c)). The dependence of this distance on the compression parameters

*G*and Δ

*x*is also shown in Fig. 23(d). There exists a strong correlation between the accuracy of the FMM computations and the distance between physical and numerical surface. As a matter of fact, the results of the FMM computations for the nondifferentiable mesh are optimal when the physical and the numerical surface coincide and deteriorate rapidly, when we move away from these sweet spots.

*τ*= 0.002. As expected, the above effect is visible, too, i.e., the error oscillates as a function of Δ

*x*for a fixed value of

*G*. However, and in particular for the low-lying modes, the amplitude of the oscillations is significantly reduced relative to the amplitude of the nondifferentiale mesh (cf. Fig. 23(b)); for our low-index fiber system about one order of magnitude.

*G*larger and smaller than those shown in Fig. 23, for a larger and smaller number of plane waves, and for larger values of the dielectric permittivity contrast. In all cases, the behavior as depicted in Fig. 23 is qualitatively reproduced. This leads us to our final conclusion regarding ASR and AC within FMM for dielectric structures. The nondifferentiable mesh performs well for small values of

*G*as this corresponds to a good representation of the singularities in the effective permittivity due to an increased density of coordinate lines at the surface. However, special care has to be exercised when choosing the real space discretization and the parameters for the compression function(s). More precisely, the number of real space points and the parameters

*G*and Δ

*x*have to be chosen such as to best represent the surface of the physical system. Here, an inconsistent choice of the mesh parameters easily leads to a rather significant loss of accuracy and/or erratic behavior with regard to convergence. The differentiable mesh exhibits the same qualitative behavior but is considerably less sensitive to the actual parameter choice, especially in terms of accuracy— good performance is achieved over a wide range of compression parameters

*G*and Δ

*x*. Overall, the nondifferentiable mesh with optimal parameters does deliver somewhat better results than the differentiable mesh. This is the manifestation of the grid-aligned effective permittivity of the nondifferentiable mesh which is more compliant with the Fourier factorization rules than the differentiable mesh (see section 5).

### 9.2. Metallic structures

**17**, 8051–8061 (2009). [CrossRef] [PubMed]

*d*= 700 nm) of metallic cylinders (height 50 nm, radius 150 nm) in air. The cylinders are centered in the unit cell and their axes are oriented in propagation direction. For the metal we use the permittivity given by a Drude model with the parameters

*ε*

_{∞}= 9.0685, a plasma frequency

*ω*= 1.3544 · 10

_{D}^{16}Hz and a damping coefficient

*γ*= 1.1536 · 10

^{14}Hz, corresponding to gold [15

15. A. Vial, A.-S. Grimault, D. Macías, D. Barchiesi, and M. Lamy de la Chapelle, “Improved analytical fit of gold dispersion: Application to the modeling of extinction spectra with a finite-difference time-domain method,” Phys. Rev. B **71**, 085416 (2005). [CrossRef]

*G*= 0.02 and Δ

*x*= 0.5. We find the resonance frequency at 829 nm which is the subject of our further investigation. Specifically, in Fig. 24(b) we compare the convergence behavior of the reflectance for the nondifferentiable mesh, two smoothed, and two differentiable meshes using the same compression function as in Fig. 24(a).

*G*and Δ

*x*. Consequently, we display in Fig. 25 the dependence of the reflectance at the resonance on

*G*and Δ

*x*for a fixed number of 997 plane waves and a real space discretization of 1024 × 1024 points when using different meshes. For the nondifferentiable mesh, we observe the same characteristics as in the dielectric case (smooth dependence on

*G*with a preference for low values, oscillatory behavior as function of Δ

*x*). In contrast, the differentiable mesh exhibits a rather erratic behavior over much of the parameter space. Consequently, for metallic structures we obtain (not unexpectedly) that an accurate representation of the surface is of paramount importance and that this can be facilitated best via the nondifferentiable mesh.

## 10. Conclusion

## Acknowledgments

## References and links

1. | K. Busch, G. von Freymann, S. Linden, S. F. Mingaleev, L. Tkeshelashvili, and M. Wegener, “Periodic nanostructures for photonics,” Phys. Rep. |

2. | L. Li, “Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings,” J. Opt. Soc. Am. A |

3. | P. Lalanne and G. M. Morris, “Highly improved convergence of the coupled-wave method for TM polarization,” J. Opt. Soc. Am. A |

4. | G. Granet and B. Guizal, “Efficient implementation of the coupled-wave method for metallic lamellar gratings in TM polarization,” J. Opt. Soc. Am. A |

5. | L. Li, “Use of Fourier series in the analysis of discontinuous periodic structures” J. Opt. Soc. Am. A |

6. | T. Weiss, G. Granet, N. A. Gippius, S. G. Tikhodeev, and H. Giessen, “Matched coordinates and adaptive spatial resolution in the Fourier modal method,” Opt. Express |

7. | G. Granet, “Reformulation of the lamellar grating problem through the concept of adaptive spatial resolution,” J. Opt. Soc. Am. A |

8. | G. Granet and J.-P. Plumey, “Parametric formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. A |

9. | T. Vallius and M. Honkanen, “Reformulation of the Fourier modal method with adaptive spatial resolution: application to multilevel profiles,” Opt. Express |

10. | S. Essig and K. Busch, “Generation of adaptive coordinates and their use in the Fourier Modal Method,” Opt. Express |

11. | L. Li, “Fourier modal method for crossed anisotropic gratings with arbitrary permittivity and permeability tensors,” J. Opt. A |

12. | L. Li, “New formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. Soc. Am. A |

13. | 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 |

14. | A. W. Snyder and J. D. Love |

15. | A. Vial, A.-S. Grimault, D. Macías, D. Barchiesi, and M. Lamy de la Chapelle, “Improved analytical fit of gold dispersion: Application to the modeling of extinction spectra with a finite-difference time-domain method,” Phys. Rev. B |

**OCIS Codes**

(050.1970) Diffraction and gratings : Diffractive optics

(050.1755) Diffraction and gratings : Computational electromagnetic methods

(160.3918) Materials : Metamaterials

(160.5298) Materials : Photonic crystals

**ToC Category:**

Diffraction and Gratings

**History**

Original Manuscript: May 10, 2012

Revised Manuscript: June 25, 2012

Manuscript Accepted: July 3, 2012

Published: July 16, 2012

**Citation**

Jens Küchenmeister, Thomas Zebrowski, and Kurt Busch, "A construction guide to analytically generated meshes for the Fourier Modal Method," Opt. Express **20**, 17319-17347 (2012)

http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-20-16-17319

Sort: Year | Journal | Reset

### References

- K. Busch, G. von Freymann, S. Linden, S. F. Mingaleev, L. Tkeshelashvili, and M. Wegener, “Periodic nanostructures for photonics,” Phys. Rep.444, 101–202 (2007). [CrossRef]
- L. Li, “Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings,” J. Opt. Soc. Am. A13, 1024–1034 (1996). [CrossRef]
- P. Lalanne and G. M. Morris, “Highly improved convergence of the coupled-wave method for TM polarization,” J. Opt. Soc. Am. A13, 779–784 (1996). [CrossRef]
- G. Granet and B. Guizal, “Efficient implementation of the coupled-wave method for metallic lamellar gratings in TM polarization,” J. Opt. Soc. Am. A13, 1019–1023 (1996). [CrossRef]
- L. Li, “Use of Fourier series in the analysis of discontinuous periodic structures” J. Opt. Soc. Am. A13, 1870–1876 (1996). [CrossRef]
- T. Weiss, G. Granet, N. A. Gippius, S. G. Tikhodeev, and H. Giessen, “Matched coordinates and adaptive spatial resolution in the Fourier modal method,” Opt. Express17, 8051–8061 (2009). [CrossRef] [PubMed]
- G. Granet, “Reformulation of the lamellar grating problem through the concept of adaptive spatial resolution,” J. Opt. Soc. Am. A16, 2510–2516 (1999). [CrossRef]
- G. Granet and J.-P. Plumey, “Parametric formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. A4, S145–S149 (2002). [CrossRef]
- T. Vallius and M. Honkanen, “Reformulation of the Fourier modal method with adaptive spatial resolution: application to multilevel profiles,” Opt. Express10, 24–34 (2002). [PubMed]
- S. Essig and K. Busch, “Generation of adaptive coordinates and their use in the Fourier Modal Method,” Opt. Express18, 23258–23274 (2010). [CrossRef] [PubMed]
- L. Li, “Fourier modal method for crossed anisotropic gratings with arbitrary permittivity and permeability tensors,” J. Opt. A5, 345–355 (2003). [CrossRef]
- L. Li, “New formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. Soc. Am. A14, 2758–2767 (1997). [CrossRef]
- 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. A24, 2880–2890 (2007). [CrossRef]
- A. W. Snyder and J. D. LoveOptical Waveguide Theory, (Chapman and Hall, 1983).
- A. Vial, A.-S. Grimault, D. Macías, D. Barchiesi, and M. Lamy de la Chapelle, “Improved analytical fit of gold dispersion: Application to the modeling of extinction spectra with a finite-difference time-domain method,” Phys. Rev. B71, 085416 (2005). [CrossRef]

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

### Figures

Fig. 1 |
Fig. 2 |
Fig. 3 |

Fig. 4 |
Fig. 5 |
Fig. 6 |

Fig. 7 |
Fig. 8 |
Fig. 9 |

Fig. 10 |
Fig. 11 |
Fig. 12 |

Fig. 13 |
Fig. 14 |
Fig. 15 |

Fig. 16 |
Fig. 17 |
Fig. 18 |

Fig. 19 |
Fig. 20 |
Fig. 21 |

Fig. 22 |
Fig. 23 |
Fig. 24 |

Fig. 25 |
||

« Previous Article | Next Article »

OSA is a member of CrossRef.