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

doi:10.1364/OE.20.017319

« Show journal navigation

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

Jens Küchenmeister, Thomas Zebrowski, and Kurt Busch »View Author Affiliations

Optics Express, Vol. 20, Issue 16, pp. 17319-17347 (2012)
http://dx.doi.org/10.1364/OE.20.017319

View Full Text Article

Acrobat PDF (12084 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Article Navigation

▸ Article Outline
– Abstract
1. Introduction
2. Covariant formulation of the Fourier Modal Method with generalized coordinates
3. Construction of nondifferentiable meshes
4. Construction of smoothed meshes
5. Construction of differentiable meshes
6. Enforcing periodicity
7. Nonrectangular unit cells
8. Compression of coordinate lines - adaptive spatial resolution
9. Suitable parameter choice and convergence characteristics
10. Conclusion
– References and links

Article Tools
Full Text
Article Info
References
Cited By
Figures
Metrics
Download PDF

Viewing Equations in the
HTML Article

Optimal Viewing Recommendations

Full Text
Article Info
References (15)
Cited By
Figures (25)
Metrics

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.

1. Introduction

In recent years, periodic nanostructures have gathered a tremendous amount of interest [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]

]. Both photonic crystals and metamaterials exhibit remarkable optical properties—complete photonic bandgaps, left-handed behavior, and high polarization rotation represent prominent examples. The experimental realization of such system remains a challenging field of research, so that numerical modeling is necessary both for a deeper understanding as well as for finding optimized and/or robust designs. Especially the latter requires fast and versatile numerical methods in order to cover large parameter spaces.

One of the methods used to predict the transmission properties of periodic photonic systems (both dielectric and metallic) is the Fourier Modal Method (FMM). Within FMM, one usually considers systems that are finite in 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

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]

A significant advancement has been achieved when Lalanne and Morris [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]

] and Granet and Guizal [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]

] found an improved formulation of the FMM for lamellar gratings. In fact, it is crucial for the FMM that the permittivity distribution is correctly Fourier transformed. The corresponding Fourier factorization rules have been put on a rigorous mathematical footing by Li [5

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

]. However, calculating the response of a metallic structure remained difficult due to in-plane stair-casing effects for not-grid-aligned permittivity distributions and most of all the Gibbs phenomenon. The latter refers to the oscillations that occur when a step function such as the permittivity distribution is approximated by a finite Fourier series.

Two interrelated approaches are pursued to address these challenges: Adaptive Coordinates (AC) and Adaptive Spatial Resolution (ASR). Both concepts employ coordinate transformations that lead to modifications of the permittivity distribution so that the convergence of FMM is accelerated. Within the AC approach, the coordinate lines are adapted to the surface of the structure [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]

]. Thereby, the structures exhibit straight and grid-aligned surfaces in the transformed space. The coordinate transformations for ASR increase the density of coordinate lines close to the structure surface so that discontinuities in the permittivity are resolved better and the problems associated with the Gibbs phenomenon are alleviated [7

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]

, 8

G. Granet and J.-P. Plumey, “Parametric formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. A 4, S145–S149 (2002). [CrossRef]

, 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]

Hence, the challenge is to find suitable coordinate transformations. During the past few years, two conceptually different approaches have emerged: in Ref. [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]

] an automated, numerical generation of meshes was presented. This offers a versatile tool since arbitrary permittivity distributions can be treated. The second approach is to find analytical coordinate transformations such as in Ref. [6

]. Here, an analytical coordinate transformation for a square array of cylinders has been presented. Such analytical mesh generation offers a lot of freedom in terms of selective and systematic change in the mesh properties. It can also be very time-efficient: when a large parameter scan is required in which the geometrical parameters of the structure are varied, analytical transformations can be easily adjusted.

In this work, we focus on how analytical meshes can be constructed even for complicated structures. In addition, we investigate the role of differentiability of the meshes and which mesh parameters should be used for a specific structure.

In section 2, we briefly recapitulate the basics of the FMM approach with AC and ASR. We follow up in sections 3, 4 and 5 with descriptions how to construct non-differentiable, smoothed, and differentiable AC meshes, respectively. In sections 6 and 7, we provide details on the realization of periodic meshes for complex geometries and non-rectangular unit cells. Subsequently, in section 8, we describe the compression of coordinate lines of AC meshes in order to realize ASR. Finally, in section 9, we carry out a comprehensive study regarding the parameter choices and the convergence characteristics for the different meshes. Based on this, we formulate guidelines for the choice of the mesh and corresponding parameters for specific structures.

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

In this section, we discuss the basics of the FMM in generalized coordinates. A typical system of interest is displayed in Fig. 1.

Fig. 1 Schematic view of a typical system of interest in FMM computations. The system is finite in x̄³-direction and periodic with respect to the x̄¹x̄² plane. k_in denotes the wave vector of an incident plane wave. Ox̄¹x̄²x̄³ describes a Cartesian coordinate system.

In content and notation we follow previous publications on the topic [6

, 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]

] so that we may be brief. Explicitly, for a given layer we distinguish between a Cartesian coordinate system Ox̄¹x̄²x̄³ and a curvilinear coordinate system Ox¹x²x³. They are connected by a transformation within the lateral plane of the general form

{\bar{x}}^{1} = {\bar{x}}^{1} (x^{1}, x^{2}),

(1)

{\bar{x}}^{2} = {\bar{x}}^{2} (x^{1}, x^{2}),

(2)

{\bar{x}}^{3} = x^{3} .

(3)

Our aim is to solve Maxwell’s curl equations which are given in covariant form as [11

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

]

ξ^{ρ α τ} \partial_{σ} E_{τ} = i k_{0} \sqrt{g} μ^{ρ σ} H_{σ},

(4)

ξ^{ρ σ τ} \partial_{σ} H_{τ} = - i k_{0} \sqrt{g} ε^{ρ σ} E_{σ} .

(5)

In these equations, ξ 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₀ = ω/c is the vacuum wave number (c denotes the vacuum speed of light and ω the frequency). The metric tensor is

g^{ρ σ} = \frac{\partial x^{ρ}}{\partial {\bar{x}}^{τ}} \frac{\partial x^{σ}}{\partial {\bar{x}}^{τ}},

(6)

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

ε^{ρ σ} = \frac{\partial x^{ρ}}{\partial {\bar{x}}^{τ}} \frac{\partial x^{σ}}{\partial {\bar{x}}^{χ}} {\bar{ε}}^{τ χ},

(7)

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.

From here onward, we obtain a working scheme by proceeding along the lines of standard FMM—we employ a Fourier-Floquet ansatz for the fields, including a plane wave ansatz in x̄³ direction, and formulate an eigenvalue problem using the correct Fourier factorization rules [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]

]. The fields are expanded into the corresponding eigenmodes and the procedure is carried out for all layers. The layers are connected using a scattering matrix which ensures boundary conditions between adjacent layers [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]

]. This gives the expansion coefficients for the fields in each layer. Finally, we transform the resulting fields in the corresponding regions back into the Cartesian system in order to find the transmittance and reflectance coefficients [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]

]. Since we focus in this manuscript on the construction of the transformations (1) and (2), we refer to Ref. [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]

] for further details on the FMM with and without coordinate transformation. For an alternative approach based on normal-to-the-surface vectors, see Ref. [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]

Before we present in the subsequent sections different approaches to the actual construction of coordinate transformations, we describe the general idea of adaptive coordinates in more detail by way of an example. Without coordinate transformation, the permittivity distribution within a unit cell of a cylindrical fiber array may look like the one depicted in Fig. 2(a). Here, we display the material matrix of a fiber with ε_fib = 2 in a background material with ε_bg = 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

] and a corresponding sector of the resulting mesh is shown in Fig. 2(b). The mesh has two important properties. First, there are many points where the mesh is not differentiable—some of them are marked with green, dashed circles. Second, there is one point (marked with a red circle) where the coordinate lines in x̄¹- and in x̄²-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

\sqrt{g} ε^{ρ σ}

, with the square root of the reciprocal of the metric determinant

\sqrt{g}

and ε^ρσ like in Eq. (7). Therefore, we refer to

\sqrt{g} ε^{ρ σ}

from now on as the effective permittivity. For our present example, we depict the

\sqrt{g} ε^{11}

-component of the tensor in Fig. 2(c). Several features are worth pointing out: First, the structure is now grid-aligned, i.e., the coordinate transformation maps the circular fiber cross section onto a squarish geometry for which the Fourier factorization rules are much more effective. Second, depending on the tensor component under investigation the nondifferentiable points in the mesh lead to discontinuities in the effective permittivity. The third feature is actually not visible here since the color scale is saturated at 4: In the corners of the material square, values in excess of 5000 are reached. In order to understand why, we investigate the derivative ∂x¹/∂x̄¹. It can be obtained by (for other derivatives, the situation is completely analogous)

\frac{\partial x^{1}}{\partial {\bar{x}}^{1}} = \underset{= {| J |}^{- 1}}{\underset{︸}{{(\frac{\partial {\bar{x}}^{1}}{\partial x^{1}} \frac{\partial {\bar{x}}^{2}}{\partial x^{2}} - \frac{\partial {\bar{x}}^{1}}{\partial x^{2}} \frac{\partial {\bar{x}}^{2}}{\partial x^{1}})}^{- 1}}} \frac{\partial {\bar{x}}^{2}}{\partial x^{2}} .

(8)

If the coordinate lines of the mesh are parallel in x̄¹- and x̄²-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¹/∂x̄¹ and, thereby, to a diverging effective permittivity.

Fig. 2 Illustration of the effect of adaptive coordinates. Panel (a) depicts the permittivity distribution for a fiber with radius 0.3 centered in a square unit cell within a Cartesian coordinate system. Panel (b) displays a section of the nondifferentiable fiber mesh from Ref. [6

]. Some points where the mesh is nondifferentiable have been marked with green, dashed circles. The point where the coordinate lines of the mesh are parallel in x̄¹- and x̄²-direction is marked with a red circle. Panel (c) displays the

\sqrt{g} ε^{11}

-component of the effective permittivity tensor that is obtained when the adaptive coordinates from (b) are applied to (a). In the transformed space, the effective permittivity tensor is fully grid-aligned but the numerical values at the rectangle edges are several orders of magnitude larger than 4. The color scale has been saturated at 4 in order to show more features.

In the following sections, we show how to construct various classes of meshes. The first class is very similar to the one just discussed—in the transformed space, the material distribution assumes a rectangular form. As hinted to above, this means that the material surface is enclosed by four specific coordinate lines and no specific attention to smoothness is paid. We delineate how to construct such 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.

At the end of this section, we would like to point out the two distinct interpretations that may be associated with the coordinate transformations discussed here. The first interpretation is that we map a given permittivity distribution onto a new distribution according to a given mesh. This means that we go from bent coordinate lines to straight coordinate lines (or material surfaces) as shown in Fig. 2. This interpretation is particularly useful when we discuss the shape of the transformed material matrix in the FMM. The second interpretation is used for constructing the meshes. Then, we talk about how to map straight coordinate lines onto bent lines that match the surface given by the material distribution. Both these interpretations are valid—which one to use depends on whether the construction or the application of the meshes is emphasized.

3. Construction of nondifferentiable meshes

In this section, we show how to construct nondifferentiable meshes associated with complex structures by way of a particular example, the crescent depicted in Fig. 3(a). Its tips are rounded so that it is defined by six quantities: the radius r₁ and the center (M₁, N₁) of the large (outer) circle, the radius r₂ and the center (M₂, N₂) of the small (inner) circle, and the radii r₃ and r₄ 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₄ = r₃) and not rotated in the unit cell, i.e., M₁ = M₂. 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.

Fig. 3 Layout of a complex unit cell for the illustration of mesh construction. A crescent with rounded tips is described by four circles—a large, outer circle (blue, radius r₁), a small, inner circle (red, radius r₂), and two tiny circles at the crescent tips (radii r₃ and r₄). To define the orientation of the crescent within the unit cell, the centers (M₁, N₁) and (M₂, N₂) of the inner and outer circles need to be specified. The centers (M₃, N₃) and (M₄, N₄) of the tiny circles are uniquely determined, once the centers of the inner and outer circles and the radii of all four circles are specified.

The crescent is composed of four circle arcs described by

{(x^{1} - M_{i})}^{2} + {(x^{2} - N_{i})}^{2} = r_{i}^{2}, i = 1, \dots, 4,

(9)

where the circles are labeled by their numbers 1 to 4. After some geometric considerations (see Fig. 3(b)), we find the center of circle 3 at

N_{3} = \frac{{(r_{2} + r_{3})}^{2} - {(r_{1} - r_{3})}^{2} + N_{1}^{2} - N_{2}^{2}}{2 N_{1} - 2 N_{2}}, M_{3} = M_{1} - \sqrt{{(r_{1} - r_{3})}^{2} - {(N_{3} - N_{1})}^{2}} .

(10)

The point

P = (x_{P}^{1}, x_{P}^{2})

where the tip circle 3 meets the outer circle 1 is given by

x_{P}^{1} = M_{3} - r_{3} cos (arctan (\frac{N_{3} - N_{1}}{M_{1} - M_{3}})), x_{P}^{2} = N_{3} + \sqrt{r_{3}^{2} - {(x_{P}^{1} - M_{3})}^{2}},

(11)

while the intersection point

D = (x_{D}^{1}, x_{D}^{2})

of circle 3 and 2 is

x_{D}^{1} = M_{3} + r_{3} cos (arctan (\frac{N_{2} - N_{3}}{M_{2} - M_{3}})), x_{D}^{2} = N_{3} + \sqrt{r_{3}^{2} - {(x_{D}^{1} - M_{3})}^{2}} .

(12)

We now proceed with the actual meshing procedure. First, as many coordinate lines as necessary are selected so that their mapping to the structure’s surface leads to its complete coverage. Second, the surface is parametrized so as to obtain an analytical expression for the mapping of the coordinate lines from step one. Third, the mappings for the remaining coordinate lines are obtained by linear interpolation between the selected coordinate lines.

For the linear transition we define the useful function

L T (c, \bar{c}, d, \bar{d}, x) = \frac{\bar{d} - \bar{c}}{d - c} x + \bar{c} - c \frac{\bar{d} - \bar{c}}{d - c}

(13)

which, as a function of 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).

Fig. 4 Illustration of a linear transition from Cartesian coordinate lines to curved coordinate lines via the linear transition function LT. The straight coordinate lines a and b are, respectively, mapped onto the green and blue bent coordinate lines. The mapped lines in between (marked with crosses) are given by the linear transition between the bent outer lines.

Based on this, we illustrate in Figs. 5(a) and 5(b) the meshing of the crescent-shaped structure introduced in Fig. 3. First, the four characteristic points P, Q, R and S are found. Point P is given by Eq. (11) and point

R = (x_{R}^{1}, x_{R}^{2})

is given by the intersection of the (Cartesian) coordinate line

x_{P}^{1}

with the structure surface. Then, points Q and R follow from the structure’s symmetry.

R = (x_{P}^{1}, N_{1} - \sqrt{r_{1}^{2} - {(x_{R}^{1} - M_{1})}^{2}}), Q = (2 M_{1} - x_{P}^{1}, x_{P}^{2}), S = (x_{Q}^{1}, x_{R}^{2}) .

(14)

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̄¹(x¹, x²), x̄² (x¹, x²)) 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̄² (corresponding to the horizontal red lines) is the identity mapping

{\bar{x}}^{2} (x^{1}, x^{2}) = x^{2}, (x^{1}, x^{2}) \in ④ .

(15)

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

\begin{array}{l} {\bar{x}}^{1} (x^{1}, x^{2}) & = & L T (\underset{map left edge onto itself}{\underset{︸}{0, 0}}, & \underset{map P R coordinate line on circle arc}{\underset{︸}{x_{P}^{1}, M_{1} - \sqrt{r_{1}^{2} - {(x^{2} - N_{1})}^{2}}}}, x^{1}) \\ \underset{=}{(13)} & M_{1} - \sqrt{r_{1}^{2} - {(x^{2} - N_{1})}^{2}} \frac{x^{1}}{x_{P}^{1}}, & (x^{1}, x^{2}) \in ④ . \end{array}

(16)

The interior coordinate lines have thus been mapped with the help of the linear transition LT as described above. Analogous mappings can be constructed for all other zones in the unit cell [0, 1] × [0, L].

Fig. 5 Illustration of how to mesh complex structures. Specific coordinate lines that pass through the structure’s characteristic points are selected (panel (a)) and subsequently mapped onto the structure’s surface.

A compact representation of the complete coordinate transformation is facilitated by defining circle arc (CA) functions for either the left, right, top, or bottom side of the crescent. Their form is derived from Eq. (9) by choosing the correct sign of the square root

\begin{matrix} \begin{matrix} C A_{L / R} (i, x^{2}) = M_{i} \mp \sqrt{r_{i}^{2} - {(x^{2} - N_{i})}^{2}}, \\ C A_{T / B} (i, x^{1}) = N_{i} \pm \sqrt{r_{i}^{2} - {(x^{1} - M_{i})}^{2}}, \end{matrix} & i = 1, \dots, 4. \end{matrix}

(17)

Here, the subscripts L, R, T, and B refer to left, right, top, and bottom, respectively.

Finally, upon introducing the point

E = (x_{E}^{1}, x_{E}^{2}) = (2 M_{1} - x_{D}^{1}, x_{D}^{2})

, we arrive at the complete coordinate transformation as

{\bar{x}}^{1} (x^{1}, x^{2}) = {\begin{array}{l} x^{1} & (x^{1}, x^{2}) \in ①, ②, ③, ⑦, ⑧, ⑨ \\ L T (0, 0, x_{P}^{1}, C A_{L} (1, x^{2}), x^{1}) & (x^{1}, x^{2}) \in ④ \\ L T (x_{P}^{1}, C A_{L} (1, x^{2}), x_{Q}^{1}, C A_{R} (1, x^{2}), x^{1}) & (x^{1}, x^{2}) \in ⑤ \\ L T (x_{Q}^{1}, C A_{R} (1, x^{2}), 1, 1, x^{1}) & (x^{1}, x^{2}) \in ⑥ . \end{array}

(18)

The x̄² mapping is more complex but the construction follows the same principles. It reads

\begin{matrix} {\bar{x}}^{2} (x^{1}, x^{2}) = \\ {\begin{matrix} x^{2} & (x^{1}, x^{2}) \in ①, ③ ​, ④, ⑥, ⑦, ⑨ \\ \begin{matrix} L T (L, L, x_{P}^{2}, C A_{T} (3, x^{1}), x^{2}) & x^{1} \leq x_{D}^{1} \\ L T (L, L, x_{P}^{2}, C A_{B} (2, x^{1}), x^{2}) & x^{1} \in (x_{D}^{1}, x_{E}^{1}) \\ L T (L, L, x_{P}^{2}, C A_{T} (4, x^{1}), x^{2}) & x_{E}^{1} \leq x^{1} \end{matrix}} & (x^{1}, x^{2}) \in ② \\ \begin{matrix} L T (x_{P}^{2}, C A_{T} (3, x^{1}), x_{R}^{2}, C A_{B} (1, x^{1}), x^{2}) & x^{1} \leq x_{D}^{1} \\ L T (x_{P}^{2}, C A_{B} (2, x^{1}), x_{R}^{2}, C A_{B} (1, x^{2}), x^{2}) & x^{1} \in (x_{D}^{1}, x_{E}^{1}) \\ L T (x_{P}^{2}, C A_{T} (4, x^{1}), x_{R}^{2}, C A_{B} (1, x^{1}), x^{2}) & x_{E}^{1} \leq x^{1} \end{matrix}} & (x^{1}, x^{2}) \in ⑤ \\ L T (x_{R}^{2}, C A_{B} (1, x^{1}), 0, 0, x^{2}) & (x^{1}, x^{2}) \in ⑧ . \end{matrix} \end{matrix}

(19)

In Fig. 6(a), we display the resulting mesh. In addition, we may compress the coordinate lines towards the crescent’s surface as depicted in Fig. 6(b). However, we postpone a discussion of coordinate line compression to section 8.

Fig. 6 Nondifferentiable meshes for a crescent with parameters r₁ = 0.25, r₂ = 0.2, r₃ = r₄ = 0.02, M₁ = 0.45, N₁ = 0.38, M₂ = M₁, N₂ = 0.5884, L = 0.8. Panel (a) displays the basic analytical mesh. Panel (b) depicts a mesh where an additional coordinate line compression using the parameters G = 0.05, Δx¹ = 0.4,

{\bar{x}}_{1}^{1} = x_{P}^{1}

{\bar{x}}_{2}^{1} = x_{Q}^{1}

, Δx² = 0.4,

{\bar{x}}_{1}^{2} = x_{R}^{2}

{\bar{x}}_{2}^{2} = x_{P}^{2}

has been applied (see section 8 for details). Panel (c) features a close-up of (a) of the left crescent tip.

The construction principle described above represents a very general procedure and can be applied to a great variety of structures. For instance, in Fig. 7 we display characteristic points, selected coordinate lines and the resulting mesh for a step-index fiber.

Fig. 7 Characteristic points, selected coordinate lines, and nondifferentiable mesh for a step-index fiber. The outer radius is 0.3 and the inner radius is 0.1. The circles are centered in (0.45,0.4) and the unit cell is [0, 1] × [0, 0.8].

A similar mesh for a circular structure has been developed by Weiss et al. [6

] (see also Fig. 2(b)). In the subsequent sections, we will use this simple structure in order to illustrate the construction of smoothed and fully differentiable meshes and to compare their performance in actual FMM computations.

4. Construction of smoothed meshes

The construction principle for the above meshes has two key properties. First, at the characteristic points the meshes are typically nondifferentiable. Second, at some points the coordinate lines are parallel in x̄¹- and in x̄²-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.

For illustrative purposes, we consider a square unit cell [0, 1] × [0, 1] with a circle of radius 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

x_{\pm} = 0.5 \pm \frac{r}{\sqrt{2}}

In order to obtain partial differentiability for the characteristic coordinate lines, we smooth the transition from straight line to circle arc with a parabola (see Fig. 8(a)). The parabola is adjusted such that it is differentiable both at point 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

g (x^{1}, a, b, x_{-}) = b {(x^{1} - a)}^{2} + x_{-}

(20)

for the parabola g we assure continuity and differentiability at a. Demanding continuity and differentiability at ū determines the parameters a and b.

a = - \frac{2}{\bar{u} - 0.5} (- \sqrt{r^{2} - {(\bar{u} - 0.5)}^{2}} + 0.5 - x_{-}) \sqrt{r^{2} - {(\bar{u} - 0.5)}^{2}} + \bar{u},

(21)

b = \frac{1}{4} {(\bar{u} - 0.5)}^{2} {[(- \sqrt{r^{2} - {(\bar{u} - 0.5)}^{2}} + 0.5 - x_{-}) (r^{2} - {(\bar{u} - 0.5)}^{2})]}^{- 1} .

(22)

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

C A_{\pm} (x^{1}) = \pm \sqrt{r^{2} - {(x^{1} - 0.5)}^{2}} + 0.5

. The mapping for the horizontal lines, i.e., the component x̄² then reads

\begin{matrix} {\bar{x}}^{2} (x^{1}, x^{2}) = \\ {\begin{array}{l} \begin{matrix} x^{2}, & x^{1} \in [0, a) \cup (1 - a, 1], \end{matrix} & x^{2} \in [0, 1] \\ \begin{matrix} L T (0, 0, x_{-}, g (x^{1}, a, b, x_{-}), x^{2}), & x^{1} \in [a, \bar{u}] \\ L T (0, 0, x_{-}, C A_{-} (x^{1}), x^{2}) & x^{1} \in (\bar{u}, 1 - \bar{u}) \\ L T (0, 0, x_{-}, g (x^{1}, 1 - a, b, x_{-}), x^{2}), & x^{1} \in [1 - \bar{u}, 1 - a] \end{matrix}} & x^{2} \in [0, x_{-}] \\ \begin{matrix} L T (x_{-}, g (x^{1}, a, b, x_{-}), x_{+}, g (x^{1}, a, - b, x_{+}), x^{2}), & x^{1} \in [a, \bar{u}] \\ L T (x_{-}, C A_{-} (x^{1}), x_{+}, C A_{+} (x^{1}), x^{2}) & x^{1} \in (\bar{u}, 1 - \bar{u}) \\ L T (x_{-}, g (x^{1}, 1 - a, b, x_{-}), x_{+}, g (x^{1}, 1 - a, - b, x_{+}), x^{2}), & x^{1} \in [1 - \bar{u}, 1 - a] \end{matrix}} & x^{2} \in (x_{-}, x_{+}) \\ \begin{matrix} L T (x_{+}, g (x^{1}, a, - b, x_{+}), 1, 1, x^{2}), & x^{1} \in [a, \bar{u}] \\ L T (x_{+}, C A_{+} (x^{1}), 1, 1, x^{2}), & x^{1} \in (\bar{u}, 1 - \bar{u}) \\ L T (x_{+}, g (x^{1}, 1 - a, - b, x_{+}), 1, 1, x^{2}), & x^{1} \in [1 - \bar{u}, 1 - a] \end{matrix}} & x^{2} \in [x_{+}, 1] \end{array} \end{matrix}

(23)

Fig. 8 Illustration of the construction principle for smoothed meshes. Panel (a) depicts how the nondifferentiable transition from (Cartesian) characteristic line to a circle arc (as described in section 3) is smoothed by a (differentiable) parabola. The parabola and the circle arc meet at the point ū = x₋ + τ with smoothing parameter τ. Straight line and parabola meet at a which also depends on τ. Panels (b) and (c) show how the characteristic coordinate lines passing through P, Q, R and S are mapped.

The mapping for the vertical lines, i.e., the component x̄¹ can be found in an analogous manner. For our simple system, symmetry mandates that the coordinate x̄¹ is constructed from the horizontal line mapping, Eq. (23), by interchanging x¹ and x². 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).

Fig. 9 Smoothed meshes for a circular structure (radius r = 0.3) centered in a square unit cell. Panels (a) and (b) depict meshes with smoothing parameters τ = 0.07 and τ = 0.035, respectively. Panel (c) shows a close-up of the mesh in (b). See text for details.

Nevertheless, upon inspecting the above expressions we find that the mesh is nondifferentiable at points x¹ ∈ (a, 1 − a), x² = x_± and x¹ = x_±, x² ∈ (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

]. In the remainder, we refer to this mesh as the 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

In the case of smoothed meshes we performed the ”parabola construction” only for the coordinate lines passing through the characteristic points. The coordinate lines in between were obtained by linear interpolation and it is this linear interpolation that eventually leads to discontinuous partial derivatives. Therefore, we can obtain fully differentiable meshes by applying the ”parabola-construction” to 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.

Fig. 10 Illustration of the construction principle for differentiable meshes. Panel (a) shows the partitioning of the unit cell and panel (b) depicts the resulting differentiable mapping.

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

] or section 4); the mappings read

{\bar{x}}^{1} (x^{1}, x^{2}) = x^{1}, {\bar{x}}^{2} (x^{1}, x^{2}) = x^{2}, (x^{1}, x^{2}) \in ①,

(24)

{\bar{x}}^{1} (x^{1}, x^{2}) = \frac{x^{1}}{x_{-}} (0.5 - \sqrt{r^{2} - {(x^{2} - 0.5)}^{2}}), {\bar{x}}^{2} (x^{1}, x^{2}) = x^{2}, (x^{1}, x^{2}) \in ②

(25)

{\bar{x}}^{1} (x^{1}, x^{2}) = \frac{2 x^{1} - 1}{x_{+} - x_{-}} \sqrt{r^{2} - {(x^{2} - 0.5)}^{2}} + \frac{1}{2}, {\bar{x}}^{2} (x^{1}, x^{2}) = {\bar{x}}^{1} (x^{2}, x^{1}), (x^{1}, x^{2}) \in ③ .

(26)

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

r^{2} = {({\bar{x}}^{1} - x_{0}^{1})}^{2} δ + {({\bar{x}}^{2} - x_{0}^{2})}^{2} γ

(27)

with the parameters

x_{0}^{1}

x_{0}^{2}

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

x_{0}^{1} = \frac{1}{2} \frac{x^{1}}{x_{-}}

x_{0}^{2} = 0.5

δ = {(\frac{x_{-}}{x^{1}})}^{2}

, γ = 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

{\bar{x}}^{1} (x^{1}, x^{2}) = b (x^{1}) {(x^{2} - a)}^{2} + x^{1}, {\bar{x}}^{2} (x^{1}, x^{2}) = x^{2}, (x^{1}, x^{2}) \in ④ .

(28)

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¹) parametrizes the family of parabolas. Continuity and differentiability at x̄² = ū require that

b (x^{1}) = \frac{\frac{1}{4} {(\bar{u} = x_{0}^{2})}^{2}}{\frac{δ}{γ^{2}} (r^{2} - γ {(\bar{u} - x_{0}^{2})}^{2}) (- \sqrt{\frac{1}{δ} (r^{2} - γ {(\bar{u} - x_{0}^{2})}^{2})} + x_{0}^{1} - x^{1})} .

(29)

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̄¹(x¹, x²) component in zone ⑥ which follow directly from the transformations in the neighboring zones ② and ③. Continuity requires

{\bar{x}}^{1} (a, x^{2}) \underset{=}{(25)} \frac{a}{x_{-}} (0.5 - \sqrt{r^{2} - {(x^{2} - 0.5)}^{2}}), {\bar{x}}^{1} (\bar{u}, x^{2}) \underset{=}{(26)} \frac{2 \bar{u} - 1}{x_{+} - x_{-}} \sqrt{r^{2} - {(x^{2} - 0.5)}^{2}} + 0.5,

(30)

whereas differentiability requires

\begin{array}{l} {\frac{\partial {\bar{x}}^{1}}{\partial x^{1}} |}_{(a, x^{2})} \underset{=}{(25)} \frac{1}{x_{-}} (0.5 - \sqrt{r^{2} - {(x^{2} - 0.5)}^{2}}), & {\frac{\partial {\bar{x}}^{1}}{\partial x^{1}} |}_{(\bar{u}, x^{2})} \underset{=}{(26)} \frac{2 \sqrt{r^{2} - {(x^{2} - 0.5)}^{2}}}{x_{+} - x_{-}}, \\ {\frac{\partial {\bar{x}}^{1}}{\partial x^{2}} |}_{(a, x^{2})} \underset{=}{(25)} \frac{a}{x_{-}} \frac{x^{2} - 0.5}{\sqrt{r^{2} - {(x^{2} - 0.5)}^{2}}}, & {\frac{\partial {\bar{x}}^{1}}{\partial x^{2}} |}_{(\bar{u}, x^{2})} \underset{=}{(26)} - \frac{2 \bar{u} - 1}{x_{+} - x_{-}} \frac{x^{2} - 0.5}{\sqrt{r^{2} - {(x^{2} - 0.5)}^{2}}} . \end{array}

(31)

Next, we make an ansatz for the mapping inside zone ⑤ using a function f which is assumed to be differentiable but otherwise still unknown. The ansatz reads

{\bar{x}}^{1} (x^{1}, x^{2}) = {\bar{x}}^{1} (a, x^{2}) \frac{x^{1} - \bar{u}}{a - \bar{u}} + {\bar{x}}^{1} (\bar{u}, x^{2}) \frac{x^{1} - a}{\bar{u} - a} + f (x^{1}, x^{2}) (x^{1} - \bar{u}) (x^{1} - a) .

(32)

This ansatz fulfills the continuity requirements of Eq. (30) provided f does not exhibit a pole at x¹ = a or x¹ = ū. Further, the ansatz ensures that the derivative with respect to x² fulfills Eq. (31) provided (∂f/∂x²) does not exhibit poles at x¹ = a or x¹ = ū. The remaining requirements refer to the derivative of Eq. (32) with respect to x¹ which reads

\begin{array}{l} \frac{\partial {\bar{x}}^{1}}{\partial x^{1}} (x^{1}, x^{2}) & = & \frac{{\bar{x}}^{1} (a, x^{2})}{a - \bar{u}} + \frac{{\bar{x}}^{1} (\bar{u}, x^{2})}{\bar{u} - a} + \frac{\partial f}{\partial x^{1}} (x^{1}, x^{2}) (x^{1} - \bar{u}) (x^{1} - a) \\ + f (x^{1}, x^{2}) (x^{1} - a) + f (x^{1}, x^{2}) (x^{1} - \bar{u}) . \end{array}

(33)

Evaluating Eq. (33) at ū and a and comparing with Eq. (31) results in

\begin{array}{l} f (\bar{u}, x^{2}) = \frac{1}{\bar{u} - a} (- \frac{{\bar{x}}^{1} (a, x^{2})}{a - \bar{u}} - \frac{{\bar{x}}^{1} (\bar{u}, x^{2})}{\bar{u} - a} + {\frac{\partial {\bar{x}}^{1}}{\partial x^{1}} |}_{(\bar{u}, x^{2})}), \\ f (a, x^{2}) = \frac{1}{a - \bar{u}} (- \frac{{\bar{x}}^{1} (a, x^{2})}{a - \bar{u}} - \frac{{\bar{x}}^{1} (\bar{u}, x^{2})}{\bar{u} - a} + {\frac{\partial {\bar{x}}^{1}}{\partial x^{1}} |}_{(a, x^{2})}) . \end{array}

(34)

We would like to emphasize that every function f that fulfills Eq. (34) and does not have poles at x¹ = a or x¹ = ū 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

f (x^{1}, x^{2}) = \frac{f (\bar{u}, x^{2}) - f (a, x^{2})}{\bar{u} - a} x^{1} + f (\bar{u}, x^{2}) - \frac{f (\bar{u}, x^{2}) - f (a, x^{2})}{\bar{u} - a} \bar{u} .

(35)

Thereby, we acquired a differentiable mapping x̄¹(x¹, x²) in zone ⑤—namely Eq. (32) with f given by Eq. (35).

We obtain the remaining differentiable mapping x̄²(x¹, x²) in this zone by a ”parabola-construction” as in the case of zone ④, this time, however, in the x²-direction

{\bar{x}}^{2} (x^{1}, x^{2}) = b (x^{2}) {(x^{1} - a)}^{2} + x^{2} .

(36)

Then, b(x²) is identical to Eq. (29) except that δ is interchanged with γ and

x_{0}^{1}

with

x_{0}^{2}

. The ellipse parameters taken from the x̄² transformation in zone ③ (see Eq. (26)) are

x_{0}^{1} = 0.5, x_{0}^{2} = 0.5, γ = {(\frac{x_{+} - x_{-}}{2 x^{2} - 1})}^{2}, δ = 1.

(37)

Zone ⑥: In this zone we only need to find the mapping of one component because any point in this area obeys the symmetry relation

({\bar{x}}^{1} (x^{1}, x^{2}), {\bar{x}}^{2} (x^{1}, x^{2})) = ({\bar{x}}^{2} (x^{2}, x^{1}), {\bar{x}}^{1} (x^{2}, x^{1})) .

(38)

The boundary conditions for this zone are derived from the mappings in zones ④ and ⑤. In particular, we use symmetry to connect the x̄² mapping in zone ⑤ (Eq. (36)) to the x̄¹ mapping in the area [ū, 1−ū]×[a,ū]. With the abbreviations Δ := x₊ − x₋, v̄:= ū−0.5 and

s : = \sqrt{r^{2} - {\bar{v}}^{2}}

the boundary conditions read

\begin{array}{l} {\bar{x}}^{1} (a, x^{2}) \underset{=}{(28)} b^{'} (a) {(x^{2} - a)}^{2} + a, & with & b^{'} (a) \underset{=}{(29)} \frac{1}{4} \frac{a}{x_{-}} \frac{{\bar{v}}^{2}}{(- s + 0.5 - x_{-}) s^{2}}, \\ {\bar{x}}^{1} (\bar{u}, x^{2}) \underset{=}{(36)} b^{″} (\bar{u}) {(x^{2} - a)}^{2} + \bar{u}, & with & b^{″} (\bar{u}) \underset{=}{(37)} \frac{1}{2} \frac{{\bar{v}}^{3}}{(s - 0.5 Δ) Δ s^{2}}, \end{array}

(39)

\begin{array}{l} {\frac{\partial {\bar{x}}^{1}}{\partial x^{2}} |}_{(a, x^{2})} \underset{=}{(28)} 2 b^{'} (a) (x^{2} - a), & {\frac{\partial {\bar{x}}^{1}}{\partial x^{1}} |}_{(a, x^{2})} \underset{=}{(28)} \frac{1}{4} \frac{{\bar{v}}^{2} {(x^{2} - a)}^{2}}{(- s + 0.5 - x_{-}) x_{-} s^{2}} + 1, \\ {\frac{\partial {\bar{x}}^{2}}{\partial x^{2}} |}_{(\bar{u}, x^{2})} \underset{=}{(36)} 2 b^{″} (\bar{u}) (x^{2} - a), & {\frac{\partial {\bar{x}}^{1}}{\partial x^{1}} |}_{(\bar{u}, x^{2})} \underset{=}{(36)} \frac{1}{2} \frac{{\bar{v}}^{2} {(x^{2} - a)}^{2}}{(s - 0.5 Δ) Δ s^{2}} + 1. \end{array}

(40)

In order to find the desired mapping, we take advantage of the generality of our formalism developed above—Eqs. (32) to (35) are all perfectly valid for the present considerations. Upon inserting Eqs. (39) and (40) in Eqs. (32) to (35), we finish the construction of a mesh that is differentiable everywhere in the unit cell. The mappings in the zones are depicted separately in Fig. 10(b). The parabola transitions are recognizable very well at the edges of the mapped zones ⑤ and ⑥. The only free parameter used for the construction of the differentiable mesh of the circle is the smoothing parameter τ. In Fig. 11 we depict the resulting meshes for two different values of τ.

Fig. 11 Differentiable meshes for a circular structure (radius r = 0.3) centered in a square unit cell. Panels (a) and (b) depict meshes with smoothing parameters τ = 0.07 and τ = 0.035, respectively. Panel (c) shows a close-up of the mesh in (b). See text for details.

Finally, we find it very instructive to compare the transformed permittivities for the different mesh types. In Fig. 12, we display the

\sqrt{g} ε^{11}

-component of the effective permittivity, i.e.,

\sqrt{g} ε^{11}

as introduced in section 2, for the nondifferentiable (panel (a)), a smoothed (panel (b)), and a differentiable mesh (panel (c)) for a circular structure of radius 0.3 that is centered within a square [0, 1]×[0, 1] unit cell. The corresponding meshes are displayed in Figs. 2(b), 9(b), and 11(b), respectively. In all cases, only the area [0.2, 0.8] × [0.2, 0.8] is shown and the effective permittivity has been sampled with 1024 × 1024 points (the same sampling is used for the numerical simulations which we present in section 9). Figure 12(a) depicts the distribution we already showed in Fig. 2(c).

Fig. 12

\sqrt{g} ε^{11}

-component of the effective permittivity for different meshes for a circular structure (radius r = 0.3) that is centered in a square unit cell. Panels (a), (b), and (c) show the transformed permittivity for the nondifferentiable, the smoothed, and the differentiable meshes of Figs. 2(b), 9(b), and 11(b), respectively. For a sampling with 1024×1024 points, the values at the corners of the rectangle for the nondifferentiable mesh exceed 5000. However, we have saturated the color scale at 4 in order to be able to display more features. The smoothed and differentiable meshes have been constructed with a smoothing parameter τ = 0.035. All plots show the sector [0.2, 0.8] × [0.2, 0.8] of the [0, 1] × [0, 1] unit cell.

First, we observe that the different mappings indeed transform the circular structure to a square-shaped object. For the above sampling, the effective permittivity values for the non-differentiable mesh exceed 5000 in the corners of the square, a clear sign of the unavoidable singularities associated with applying a nondifferentiable mesh to a circular structure (see section 2). In addition, the circle is exactly mapped onto the square. For the smoothed as well as for the differentiable mesh the circle is mapped on a square with rounded corners and is, thus, not fully grid-aligned. However, the values of the effective permittivity are considerable lower in both cases. The difference between them is the slightly different behavior at the structure’s boundaries. The effective permittivity of the differentiable mesh shows jumps only at the surface of the transformed structure. In contrast, the effective permittivity of the smoothed mesh features additional discontinuities since it is not differentiable everywhere (see section 4). We find qualitatively similar behavior for all non-zero components of the effective permittivity tensor.

Based on the above observations, it is an interesting question to ask which of the meshes will exhibit the best convergence characteristics for FMM-computations. The strictly grid-aligned effective permittivity for the nondifferentiable meshes is ideally suited for the Fourier factorization rules (see section 1) but its rather high values are clearly detrimental. The situation is reversed for the smoothed and the differentiable mesh.

However, before we address this question, we first complete the discussion of the mesh construction in sections 6–8, where we will provide guidelines for the realization of periodic meshes for low-symmetry motifs in the unit cell (section 6), the construction of meshes for nonrectangular lattices (section 7), and the implementation of adaptive spatial resolution via coordinate line compression (section 8).

In order to prepare for the latter, we note in the case of the nondifferentiable and the smoothed mesh, the coordinate line that is mapped onto the circle’s surface is given by x₋ (see Fig. 8). In order to find the corresponding coordinate line α for the differentiable mesh, we numerically compute the point that fulfills

0.5 - r = {\bar{x}}^{1} (α, 0.5), (α, 0.5) \in ⑤

(41)

6. Enforcing periodicity

The meshes that we have constructed in the previous sections exhibit a mandatory requirement for being used with FMM computations: They are periodic, i.e., every coordinate line enters and leaves the unit cell at opposing unit cell edges. Apparently, certain meshes that are generated by our formalism described above do not exhibit this property (see, e.g., Fig. 13). For instance, in order to construct a mesh for a square array of ellipses, the natural choice of characteristic points are those where the major and minor axes pierce the ellipse’s surface. However, if the ellipse is rotated relative to the square lattice, our formalism leads to nonperiodic meshes.

Fig. 13 Construction points and mesh for an elliptical structure that is not aligned with the axes of the square array associated with the unit cell. This example further illustrates how our construction principle may lead to aperiodic meshes. See text for further details.

There are three ways to restore periodicity. First, one could choose a different set of characteristic points that does lead to periodic meshes. Second, one can change the way the characteristic coordinate lines are mapped. We depict this approach in Figs. 14(a) and 14(b). Here, the characteristic coordinate lines were modified by additional straight line segments (dashed) so that the coordinate lines leave the unit cell at the same point they entered at the opposite face of the unit cell. The dashed part does not have to be composed of straight lines—other functions such as higher order polynomials are conceivable as well. Third, one can use a 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.

Fig. 14 Illustration of different approaches for enforcing periodicity of meshes. Panels (a) and (b) depict the method of adding linear transitions to the outer edge of the unit cell. Panel (c) illustrates the method of the mirror structure. Color has been added to mark the actual structure. See text for further details.

In order to demonstrate that neither the mirror-structure approach nor the linear transitions to the unit cell edges lead to discernible errors, we investigate the rib waveguide sketched in Fig. 15(a). To obtain a converged result we used 2601 Fourier coefficients with a real space discretization of 1024×1024 points. The unit cell is 4μ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 ε_s = 12.1 while the background medium is air with ε_bg = 1. The incoming wavelength was chosen to be the telecom wavelength at 1550 nm.

Fig. 15 Analysis of a square array of silicon waveguides with trapezoidal cross section using the mirror-structure approach for periodic mesh construction. Panel (a) depicts a schematic of the system. Panels (b) and (c) show, respectively, the periodic mesh for the trapezoidal structure and the absolute value of the transverse electric field distribution of the fundamental mode on a logarithmic scale. There are no discernable detrimental effects due to the mirror structure that is used to ensure periodicity of the mesh. The yellow color in panel (b) and the white trapezoid in panel (c) have been added to guide the eye.

Figure 15(c) shows the absolute value of the transverse electric field of the rib waveguide’s fundamental mode. Since the permittivity is transformed according to Eq. (7), any coordinate line other than the identity will result in a distortion of the effective permittivity. Yet, Fig. 15(c) shows that the values of the numerical artifacts in the region of the mirror structure are several orders of magnitude lower than the physical fields at the trapezoid.

7. Nonrectangular unit cells

So far, we have focused on rectangular unit cells. Apparently, for many systems more freedom in the unit cell choice is highly desirable. In fact, treating nonrectangular unit cells in the framework of FMM without adaptive meshing is well known [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]

]. In this section, we will briefly discuss how adaptive meshing for a nonrectangular unit cell can be realized and present an array of circular rods in a hexagonal lattice as an example.

The basic idea is to map the hexagonal unit cell to a rectangular one (including the structure!), perform the mesh generation as discussed in the previous sections and map the rectangular unit cell (including the mesh) back to the hexagonal unit cell. A mapping between a nonrectangular unit cell with lattice angle α and a rectangular unit cell is given by [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]

]

\begin{array}{l} x_{α} = x_{rec} + y_{rec} sin α & \leftrightarrow & x_{rec} = x_{α} - y_{α} \frac{sin α}{cos α}, \\ y_{α} = y_{rec} cos α & \leftrightarrow & y_{rec} = y_{α} \frac{1}{cos α} . \end{array}

(42)

Here, the coordinates 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.

Fig. 16 Illustration of the meshing in nonrectangular unit cells. Panel (a) depicts the structure in ordinary (Cartesian) space—a circular rod in a hexagonal lattice. Panel (b) shows the transformed structure when the nonrectangular unit cell is mapped onto the rectangular unit cell—the circle turns into an ellipse. Panel (c) shows the mesh in the original unit cell that is obtained by meshing the ellipse and mapping this mesh back into ordinary (Cartesian) space. The parameters are r = 0.2, α = 30° and a [0, 1] × [0, 1] unit cell.

It is important to realize that the nonrectangular unit cell given by Eq. (42) is still described in a Cartesian coordinate system, i.e., we do not switch to an oblique-angled coordinate system. Hence, the equation of a circle remains valid and we can perform the mapping to the rectangular unit cell according to

{Circle}_{α} = {(x, y) : {(x - \frac{1}{2} (1 + sin α))}^{2} + {(y - \frac{1}{2} cos α)}^{2} = r^{2}},

(43)

\begin{array}{l} {Circle}_{rec} & \underset{=}{(42)} & {(x - y \frac{sin α}{cos α}, y \frac{1}{cos α}) : {(x - \frac{1}{2} (1 + sin α))}^{2} + {(y - \frac{1}{2} cos α)}^{2} = r^{2}}, \\ = & {(x, y) : {(x + y sin α - \frac{1}{2} (1 + sin α))}^{2} + {(y cos α - \frac{1}{2} cos α)}^{2} = r^{2}} . \end{array}

(44)

In order to perform the mesh generation analytically, we have to extract the ellipse parameters from Eq. (44). An ellipse with center (x₀, y₀) that is rotated by an angle ϕ with semi-axes c and d is described by

(\begin{matrix} x - x_{0} & y - y_{0} \end{matrix}) (\begin{matrix} cos ϕ & - sin ϕ \\ sin ϕ & cos ϕ \end{matrix}) (\begin{matrix} c^{- 2} & 0 \\ 0 & d^{- 2} \end{matrix}) (\begin{matrix} cos ϕ & sin ϕ \\ - sin ϕ & cos ϕ \end{matrix}) (\begin{matrix} x - x_{0} \\ y - y_{0} \end{matrix}) = 1.

(45)

Comparing Eq. (44) with Eq. (45) yields

c = r {(1 - sin α)}^{- 1 / 2}, d = r {(1 + sin α)}^{- 1 / 2}, ϕ = - 45 °, (x_{0}, y_{0}) = (0.5, 0.5) .

(46)

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

As stated several times above, it is often very desirable to increase the density of coordinate lines within a structure or at its surface in order to realize adaptive spatial resolution (ASR). We have used such coordinate compressions for constructing the meshes depicted in Fig. 6(b), Fig. 15(b), and Fig. 16(c). In this section, we provide the required details. In essence, ASR and adaptive coordinates (AC) are performed sequentially—first, the coordinate lines are compressed and then the mapping for the AC is applied. Technically, ASR it is just another coordinate transformation. A good proposal for a one-dimensional compression function can be found in Ref. [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]

]. In this work we have, in general, to apply two-dimensional compressions which we construct by two successive one-dimensional compressions after slightly modifying the compression functions of Ref. [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]

]. In Fig. 17(a), we display an illustrative example of such a one-dimensional compression.

Fig. 17 Illustration of how to shift the compression function from Ref. [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]

]. In all panels x is plotted on the horizontal axis and x̄ on the vertical axis. Panel (a) depicts the original compression function. Panel (b) shows how the function is shifted to larger values. The point of intersection where the transformed coordinate leaves the unit cell is denoted by x̃. In panel (c) the compression function is shifted by exactly this amount to the right. The piece that is outside the unit cell in panel (c) is attached to the left side due to periodicity. Panel (d) shows how the equidistant coordinate lines (vertical) are mapped onto compressed coordinate lines (horizontal).

The transformation function of Ref. [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]

] is

\bar{x} (x) = α + β x + \frac{γ}{2 π} sin (2 π \frac{x - x_{l - 1}}{x_{l} - x_{l - 1}}), x \in [x_{l - 1}, x_{l}], l = 2, \dots, n

(47)

with

α = \frac{x_{l} {\bar{x}}_{l - 1} - x_{l - 1} {\bar{x}}_{l}}{x_{l} - x_{l - 1}}, β = \frac{{\bar{x}}_{l} - {\bar{x}}_{l - 1}}{x_{l} - x_{l - 1}}, γ = (x_{l} - x_{l - 1}) G - ({\bar{x}}_{l} - {\bar{x}}_{l - 1}),

(48)

where {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̄_l for l = 1,..., n. Further, G denotes the slope of the function at all positions x_l. The compressed coordinates are described by x̄(x)—the transformation is a linear function superimposed with a sine. In the vicinity of the points x̄_l the coordinate line density is increased—for instance, these points could represent material interfaces along the x-coordinate line. The interval [x_l₋₁, x_l] is mapped onto the interval [x̄_l₋₁, x̄_l]. In order to maintain periodicity with period L, one chooses x̄₁ = 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

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]

] this is intended. We, on the other hand, would like to place the structures in the center of the unit cell. This amounts to shifting the compression function and we show how to accomplish this with two material interfaces. This is straightforwardly generalized to more interfaces.

We start with two parallel material interfaces in the unit cell located at two given coordinates x̄₁ and x̄₂. Furthermore, we assume that the length of the interval that is mapped onto [x̄₁, x̄₂] is given—it is referred to as Δx. Then, we perform the mapping Eq. (47) as if the first material interface lies at 0.

x \in [0, Δ x] : \bar{x} (x) = \frac{{\bar{x}}_{2} - {\bar{x}}_{1}}{Δ x} x + \frac{G Δ x - ({\bar{x}}_{2} - {\bar{x}}_{1})}{2 π} sin (2 π \frac{x}{Δ x}),

(49)

x \in [Δ x, L] : \bar{x} (x) = \tilde{α} + \tilde{β} x + \tilde{γ} sin (2 π \frac{x - Δ x}{L - Δ x}) .

(50)

Here, we have introduced the abbreviations

\tilde{α} = \frac{L ({\bar{x}}_{2} - {\bar{x}}_{1}) - Δ x L}{L - Δ x}, \tilde{β} = \frac{L - ({\bar{x}}_{2} - {\bar{x}}_{1})}{L - Δ x}, \tilde{γ} = \frac{G (L - Δ x) - (L - ({\bar{x}}_{2} - {\bar{x}}_{1}))}{2 π} .

(51)

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̄₁ upwards, see Fig. 17(b). The point where the compression function intersects with the unit cell edge, denoted x̃, is

L = \tilde{α} + \tilde{β} \tilde{x} + \tilde{γ} sin (2 π \frac{\tilde{x} - Δ x}{L - Δ x}) + {\bar{x}}_{1}

(52)

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̄₁ and x̄₂ in the interior of the unit cell then reads

x \in [0, \tilde{L}] : \bar{x} (x) = \tilde{α} + \tilde{β} (x + \tilde{x}) + \tilde{γ} sin (2 π \frac{x + \tilde{x} - Δ x}{L - Δ x}) + {\bar{x}}_{1} - L,

(53)

x \in (\tilde{L}, \tilde{L} + Δ x) : \bar{x} (x) = \frac{{\bar{x}}_{2} - {\bar{x}}_{1}}{Δ x} (x - \tilde{L}) + \frac{G Δ x - ({\bar{x}}_{2} - {\bar{x}}_{2})}{2 π} sin (2 π \frac{x - \tilde{L}}{Δ x}) + {\bar{x}}_{1},

(54)

x \in [\tilde{L} + Δ x, L] : \bar{x} (x) = \tilde{α} + \tilde{β} (x - \tilde{L}) + \tilde{γ} sin (2 π \frac{x - \tilde{L} - Δ x}{L - Δ x}) + {\bar{x}}_{1} .

(55)

We want to stress that one needs to start with the separation Δx of the interfaces instead of starting with their exact locations x₁ and x₂ (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.

In Figs. 18(a) and 18(b), we display two examples of compression functions. They are plotted together as a two dimensional mesh in Fig. 18(c) where the compression in Fig. 18(a) was applied in the horizontal direction and the compression in Fig. 18(b) in the vertical direction. When choosing the compression parameters it is important to make sure that 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̄₂ − x̄₁)/Δ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).

Fig. 18 Illustration how two one-dimensional compression functions are combined to a two-dimensional compressed mesh. The compression in (a) is applied in horizontal direction and the compression depicted in (b) (with a constant density between x̄₁ and x̄₂) is applied in the vertical direction. Panel (c) depicts the resulting mesh.

This completes our illustration of adaptive coordinates with adaptive spatial resolution: In order to obtain meshes as those depicted in Fig. 6(b), we first map a Cartesian grid to a grid that is compressed in certain regions by way of Eqs. (53) – (55). We then use this locally compressed grid as input for the adaptive coordinate mapping that deforms the grid in the desired way (nondifferentiable, smoothed, or differentiable), e.g., as in Eqs. (18) and (19).

The smoothed meshes remove the (unphysical) singularities in the effective permittivity, Eq. (7), of the nondifferentiable meshes and only certain discontinuities in the partial derivatives remain. These discontinuities, too, can be removed by using differentiable meshes. Obviously, this removal of singularities and discontinuities will be conductive for the convergence of FMM. However, another important aspect of the smoothed and the differentiable meshes is that—in contrast to the nondifferentiable meshes—the structure’s surface is no longer grid-aligned in the transformed space described by Eq. (7) (see Fig. 12). Yet, the success of the Fourier factorization rules critically depends on grid-alignment.

Which aspect is more important in which circumstances? This is the subject to be studied in the next (and final) section.

9. Suitable parameter choice and convergence characteristics

More generally, in this section we deal with the important question which type of mesh and which mesh parameters are suitable for a given structure. The number of parameters for the mesh generation is very large—for each given frequency and structure, we can vary the compression of the mesh (parameters 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.

First, we examine dielectric structures. We start with a study on how and why different parameters change the convergence results of the propagation constants of low-lying eigenmodes. This will result in guidelines how to find suitable parameters. Second, we repeat this exercise for metallic structures and their convergence characteristics.

9.1. Dielectric structures

In order to investigate the convergence behavior of the FMM with different meshes for dielectric structures, we choose a fiber with cylindrical cross section as a test system. This choice is motivated by the availability of analytical solutions for the guided eigenmodes of a single cylindrical waveguide [14

A. W. Snyder and J. D. Love Optical Waveguide Theory , (Chapman and Hall, 1983).

]. For such modes, the inherent periodicity of the FMM computations is of no relevance once the unit cell is sufficiently large, i.e., we are then performing a supercell computation. Clearly, this also means that we should exercise care for modes that are close to the cut-off as they might be fairly extended. In our subsequent computations this has been carefully checked.

Specifically, we consider a fiber of radius of 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 ε_bg = 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 λ = 800 nm.

As a measure for the mesh quality, we compute the maximum relative error of the effective refractive index of the first ten guided eigenmodes for a given number of plane waves, i.e.,

error = max {\frac{| n_{eff, analyt, 1 st mode} - n_{eff, num, 1 st mode} |}{n_{eff, analyt, 1 st mode}}, \dots, \frac{| n_{eff, analyt, 10 th mode} - n_{eff, num, 10 th mode} |}{n_{eff, analyt, 10 th mode}}},

where the effective refractive index of an eigenmode with propagation constant γ 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.

Fig. 19 Convergence characteristics of different meshes for computing the effective index of refraction of guided modes in a low-index fiber. Panel (a) displays the convergence characteristics of Cartesian, nondifferentiable, smoothed, and differentiable meshes without coordinate line compression. Panels (b) to (d) depict the corresponding dependence of the relative error (color-coded) on the compression parameters G and Δx for a fixed number of 997 plane waves. The same compression function is applied in x̄¹- and x̄²-direction. All computations have been performed on a grid with 1024 × 1024 sampling points in transformed space. Panel (b) displays the results for the nondifferentiable mesh. Panels (c) and (d) depict results for the differentiable mesh with τ = 0.002 and τ = 0.015, respectively. The sketches in panel (d) depict compression functions for four specific pairs of G and Δx. The color bar of panel (d) also applies to panels (b) and (c). See the text for further details.

As the cylindrical fiber is centered in the square unit cell we can apply the same compression function in x̄¹- and x̄²-direction. The material surface, i.e., the compression points are x̄₁ = x₋ and x̄₂ = x₊ for the nondifferentiable mesh and x̄¹ = α and x̄² = 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.

Fig. 20 Dependence of the relative error (color-coded) for a low-index fiber on the compression parameters G and Δx for 293 planes waves and 1024 × 1024 discretization points in transformed space (cf. Fig. 19 for the results of the same system with 997 plane waves). Panel (a) depicts the results for the nondifferentiable mesh and panels (b) and (c) depict the results for the differentiable mesh with τ = 0.002 and τ = 0.015, respectively. The color bar of panel (c) also applies to panels (b) and (c). See the text for further details.

Fig. 21 Dependence of the relative error (color-coded) for a high-index fiber on the compression parameters G and Δx for 997 planes waves and 1024 × 1024 discretization points (cf. Fig. 19 for the results of a low-index fiber with the same number of plane waves). Panel (a) depicts the results for the nondifferentiable mesh and panels (b) and (c) depict the results for the differentiable mesh with τ = 0.002 and τ = 0.015, respectively. The color bar of panel (c) also applies to panels (b) and (c). See the text for further details.

Fig. 22 Dependence of the relative error (color-coded) for a low-index fiber on the compression parameters G and Δx for 997 planes waves and 1000 × 1000 discretization points (cf. Fig. 19 for the results of the same system with 1024×1024 discretization points). Panel (a) depicts the results for the nondifferentiable mesh and panels (b) and (c) depict the results for the differentiable mesh with τ = 0.002 and τ = 0.015, respectively. The color bar of panel (c) also applies to panels (b) and (c). The area within the white box in panel (a) is selected for further detailed investigation in Fig. 23. See the text for further details.

From Fig. 19(b), we infer that the nondifferentiable mesh performs well for small values of 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.

To sum up our (partial) findings up to this point (see Fig. 19), nondifferentiable meshes yield the overall best results but their performance sensitively depends on the correct choice of mesh parameters. On the other hand, the differentiable meshes perform nearly as well as the optimal nondifferentiable meshes but offer the advantage of being muss less sensitive to parameter change.

As alluded to above, an interesting aspect is how the meshes perform for different numbers of plane waves. Therefore, we have performed similar computations as those depicted in Fig. 19 but this time with 293 plane waves. We display the results in Fig. 20 for the nondifferentiable mesh (panel (a)), the differentiable mesh with τ = 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.

Moreover, these findings remain valid when the dielectric contrast is increased. In Fig. 21 we display the results for computations using again 997 plane waves (as in Fig. 19) but this time for a fiber with larger permittivity ε_fib = 10. The background material is again ε_bg = 1. We have further checked that the above (partial) summary holds true for other values of the permittivity by performing calculations with varying ε_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).

An important issue which we have yet to address is the influence of the real space discretization. We have performed all the previous computations with 1024 × 1024 discretization points in order to take full advantage of the capabilities of the Fast Fourier Transform. However, we have to be aware that such a discretization must not necessarily be optimal for a given problem, including our fiber system. This may be seen from the following argument: As described in section 2, the nondifferentiable mesh leads to a diverging permittivity at specific points in space. Depending on the number of real space points we use in the computations, the coordinate lines are either close to these points or not. Hence, the maximum permittivity values ”seen” by the numerical framework changes drastically when changing the number of real space points. Another important issue is the size of the structure ”seen” by the numerics—placing a coordinate line right on the surface of the structure or next to the surface changes its effective size. In our fiber system, for instance, the effective refractive indices of the modes change according to how far away the coordinate lines are from the structure’s surface. A potential way of dealing with this problem could be increasing the number of sampling points or to find ways of appropriately distributing a given number of points. As the former option might require considerable computational resources, we will, in the following, pursue the latter option.

In order to start the quantitative analysis, we perform computations for the same low-index fiber for which we have obtained the results of Fig. 19—this time only with 1000 instead of 1024 discretization points per dimension. We depict our results in Fig. 22. As expected, the largest changes occur for the nondifferentiable mesh (cf. Fig. 19(b) and Fig. 22(a)).

Upon this (rather minimal) change of the real space discretization, the results for the differentiable meshes with τ = 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.

Indeed, while the error changes continuously with 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.

Fig. 23 Detailed investigation of the error associated with nondifferentiable and differentiable meshes. Panel (a) represents at blow-up of the results highlighted in the white box of Fig. 22(a). The maximum relative error of the first ten guided eigenmodes of a low-index fiber have been obtained with 997 plane waves and 1000 × 1000 real-space points in steps of 0.001 for G and 0.00005 for Δx. Panel (b) depicts a line-cut through (a) at G = 0.165 and the relative error of each of the first ten guided eigenmodes together with the distance of the structure’s physical surface to the numerical surface. Panel (c) illustrates the definition of the numerical surface as the middle between the two coordinate lines closest to the physical surface. Panel (d) displays the computed distance between numerical and physical surface for the same parameter range that is used in (a). Panel (e) displays the results of a similar computation as shown in (b), this time using the differentiable mesh with τ = 0.002 and G = 0.165. See text for further details.

In Fig. 23(e) we display corresponding results of the relative error of each of the first ten guided eigenmodes using a differentiable mesh with τ = 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.

We have conducted extensive studies for values of 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

We now turn to an investigation of metallic systems. While much of the above discussion qualitatively also applies in this case, the relative weighting of the above issues is rather different. This may be seen as follows. Metals are characterized by permittivities with negative real part. In turn, this leads to a finite penetration of the electromagnetic field into the metal as well as to large field enhancements near the surfaces. As a result, for metallic structures, we expect that surfaces and their adequate representation is even more important than for dielectric structures. To be specific, we consider the same system that has been discussed in Ref. [6

]. This means that we investigate a square array (lattice constant 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 ω_D = 1.3544 · 10¹⁶ Hz and a damping coefficient γ = 1.1536 · 10¹⁴ Hz, corresponding to gold [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]

]. In the absence of analytic solutions for this problem, we have computed transmittance, reflectance and absorbance spectra for normal incidence using several different compression functions and found qualitatively the same results. In Fig. 24(a) we display the corresponding spectra using the nondifferentiable mesh, 1750 plane waves and a compression with parameters 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).

Fig. 24 Convergence characteristics of different meshes for a square array of metallic cylinders of finite height. Panel (a) shows the transmittance, reflectance, and absorbance spectra that have been computed using 1750 plane waves and the nondifferentiable mesh. Our data agree well with those of Ref. [6

] for this system. Panel (b) depicts the convergence characteristics (in terms of plane waves) of the reflectance for different meshes at the resonance frequency of 829 nm. The data in panels (a) and (b) have been obtained by using the same compression function with parameters G = 0.02 and Δx = 0.5 and a real space discretization of 1024 × 1024 points.

From Fig. 24(b), we infer that the nondifferentiable mesh exhibits a considerably faster and less erratic convergence behavior than the smoothed and differentiable meshes. Along the lines of our analysis for dielectric systems, we have traced this to the fact that the transformed effective permittivity fails to be grid-aligned for the smoothed and differentiable meshes. Therefore, it is less compliant with the Fourier factorization rules.

In contrast to dielectric structures, Fig. 24(b) strongly suggests that resolving the entire surface and creating a grid-aligned structure for the transformed permittivity is crucial for the performance when investigating metals with the FMM. Here, the nondifferentiable mesh has a distinct advantage over the smoothed and the differentiable meshes.

Just as in the dielectric case, it is interesting to have a closer look at the compression parameters 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.

Fig. 25 Dependence of the on-resonance reflectance from a square array of finite-height metallic cylinders on the compression parameters G and Δx for different meshes. Panel (a) shows the results when using the nondifferentiable mesh and panel (b) shows the results when using the differentiable mesh with τ = 0.005. All computations have been performed with 997 plane waves with a real space grid of 1024×1024 sampling points. The parameter G has been stepped in 0.005 intervals and Δx in 0.0025 intervals. The color scale was saturated at 0.81. See text for further details.

10. Conclusion

In this work we have dealt with several important aspects of the FMM when using adaptive coordinates (AC) and adaptive spatial resolution (ASR). First, we have demonstrated how different issues in analytical mesh generation such as maintaining periodicity, non-rectangular unit cells and compression of coordinate lines at material interfaces can be mastered even for rather complex structures. Second, we have provided construction guidelines for different types of meshes such as nondifferentiable, smoothed and differentiable meshes. Through a careful convergence study, we have discovered ”sweet spots” for the nondifferentiable mesh where a judicious choice of real space discretization and compression parameters leads to an optimal surface representation that efficiently treats the (for this mesh unavoidable) singularities in the effective permittivity. As the nondifferentiable meshes naturally lead to the best grid-alignment, i.e., to optimal compliance with regard to the Fourier factorization rules, the nondifferentiable mesh turns out to be the best choice—provided one knows what one is doing and one can realize the correct surface representation.

The differentiable meshes exhibit a comparable performance for dielectric structures and are considerably less sensitive to the correct choice of compression parameters. This is welcome news for black-box approaches. For metals, however, obtaining grid-aligned structures is of paramount importance and the optimized nondifferentiable mesh still performs rather well (necessarily less efficient than in the dielectric case) while the smoothed and the differentiable mesh clearly fall behind (but are, of course, still much better than a simple Cartesian grid).

Acknowledgments

It is a pleasure to thank Benjamin Lutz for his work on the trapezoidal waveguide structure. In addition, we acknowledge fruitful discussions with Michael Walz and Christian Wolff. Furthermore, we thank Sabine Essig for her work on certain parts of the implementation. We acknowledge support by the Deutsche Forschungsgemeinschaft (DFG) and the State of Baden-Württemberg through the DFG-Center for Functional Nanostructures (CFN) within sub-project A1.1. T.Z. acknowledges funding by the Carl Zeiss foundation and the Karlsruhe House of Young Scientists (KHYS) and his research is embedded in the Karlsruhe School of Optics & Photonics (KSOP). Furthermore, we acknowledge support by Deutsche Forschungsgemeinschaft and Open Access Publishing Fund of the Karlsruhe Institute of Technology (KIT).

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. 444, 101–202 (2007). [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]
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]
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]
7.	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]
8.	G. Granet and J.-P. Plumey, “Parametric formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. A 4, S145–S149 (2002). [CrossRef]
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]
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]
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]
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]
14.	A. W. Snyder and J. D. Love Optical Waveguide Theory , (Chapman and Hall, 1983).
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]

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: Author | 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.

Click here to see a list of articles that cite this paper

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