The concepts of adaptive coordinates and adaptive spatial resolution significantly enhance the performance of Fourier Modal Method for the simulation of periodic photonic structures, especially metallo-dielectric systems. We present several approaches for constructing different types of analytical coordinate transformations that are applicable to a great variety of structures. In addition, we analyze these meshes with an emphasis on the resulting convergence characteristics. This allows us to formulate general guidelines for the choice of mesh type and mesh parameters.
© 2012 OSA
In recent years, periodic nanostructures have gathered a tremendous amount of interest . 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 .
A significant advancement has been achieved when Lalanne and Morris  and Granet and Guizal  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 . 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 . 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, 8, 9].
Hence, the challenge is to find suitable coordinate transformations. During the past few years, two conceptually different approaches have emerged: in Ref.  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. . 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.
In content and notation we follow previous publications on the topic [6, 10, 11] so that we may be brief. Explicitly, for a given layer we distinguish between a Cartesian coordinate system Ox̄1x̄2x̄3 and a curvilinear coordinate system Ox1x2x3. They are connected by a transformation within the lateral plane of the general form11] Eqs. (4) and (5)) denotes the reciprocal of its determinant. As a consequence of the coordinate transformation, the permittivity changes according to
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̄3 direction, and formulate an eigenvalue problem using the correct Fourier factorization rules [5, 11]. 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 . 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 . Since we focus in this manuscript on the construction of the transformations (1) and (2), we refer to Ref. [10, 11] for further details on the FMM with and without coordinate transformation. For an alternative approach based on normal-to-the-surface vectors, see Ref. .
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.  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̄1- and in x̄2-direction are parallel. In the following, we show how this transformation affects the eigenvalue problem of the FMM. The matrix that is Fourier transformed when using FMM with AC and/or ASR is not the permittivity itself but rather , with the square root of the reciprocal of the metric determinant and ερσ like in Eq. (7). Therefore, we refer to from now on as the effective permittivity. For our present example, we depict the -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 ∂x1/∂x̄1. It can be obtained by (for other derivatives, the situation is completely analogous)
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 r1 and the center (M1, N1) of the large (outer) circle, the radius r2 and the center (M2, N2) of the small (inner) circle, and the radii r3 and r4 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., r4 = r3) and not rotated in the unit cell, i.e., M1 = M2. 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.
The crescent is composed of four circle arcs described byFig. 3(b)), we find the center of circle 3 at
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 functionFig. 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).
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 is given by the intersection of the (Cartesian) coordinate line with the structure surface. Then, points Q and R follow from the structure’s symmetry.Eq. (9)). Third and final, we have to find the two-dimensional transformation function (x̄1(x1, x2), x̄2 (x1, x2)) that maps the unit cell onto itself by a linear transition between these deformed coordinate lines. This is done separately for all nine zones in the unit cell and we will discuss only zone ④ (see Figs. 5(a) and 5(b)). The red (horizontal) coordinate lines in zone ④ in Fig. 5(a) are not affected by the mapping (i.e. they are in the same position in Figs. 5(a) and 5(b)). This means that the mapping for x̄2 (corresponding to the horizontal red lines) is the identity mapping Eq. (9). Explicitly, this reads
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
Finally, upon introducing the point , we arrive at the complete coordinate transformation as
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.
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.
A similar mesh for a circular structure has been developed by Weiss et al.  (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̄1- and in x̄2-direction (see Fig. 2(b)). The latter causes a divergent effective permittivity (see Eq. (7)). Therefore, we will now describe a way to design smoothed meshes that avoid these singularities. In this context, smoothed means that the characteristic coordinate lines are mapped onto smooth curves. The result are meshes whose partial derivatives are still discontinuous at several points. In this sense, the present section represents an intermediate step from nondifferentiable to fully differentiable meshes which are discussed in the subsequent section.
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 .
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 ofFig. 8. The other coordinate lines follow by a linear transition. The functions needed are given in Eq. (13), Eq. (20) and by . The mapping for the horizontal lines, i.e., the component x̄2 then reads
The mapping for the vertical lines, i.e., the component x̄1 can be found in an analogous manner. For our simple system, symmetry mandates that the coordinate x̄1 is constructed from the horizontal line mapping, Eq. (23), by interchanging x1 and x2. 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).
Nevertheless, upon inspecting the above expressions we find that the mesh is nondifferentiable at points x1 ∈ (a, 1 − a), x2 = x± and x1 = x±, x2 ∈ (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. . 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.
Zones ①, ②, ③: In these regions, we choose the mappings exactly as in the case of the nondifferentiable circle mesh (see Ref.  or section 4); the mappings read
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 byEq. (25): , , , γ = 1. We aim to connect the straight lines (i.e., the identity mapping) in zone ① with the ellipse arcs by a family of parabolas of the form Eq. (21), b is not a constant anymore as it has been in Eq. (22)—thus b(x1) parametrizes the family of parabolas. Continuity and differentiability at x̄2 = ū require that
Zone ⑤: To find suitable differentiable mappings in the zones ⑤ and ⑥ we will utilize an alternative approach: We formulate the requirements in terms of continuity and differentiability and then make an ansatz to find appropriate functions. We start by considering the boundary conditions for the x̄1(x1, x2) component in zone ⑥ which follow directly from the transformations in the neighboring zones ② and ③. Continuity requires
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 readsEq. (30) provided f does not exhibit a pole at x1 = a or x1 = ū. Further, the ansatz ensures that the derivative with respect to x2 fulfills Eq. (31) provided (∂f/∂x2) does not exhibit poles at x1 = a or x1 = ū. The remaining requirements refer to the derivative of Eq. (32) with respect to x1 which reads Eq. (33) at ū and a and comparing with Eq. (31) results in Eq. (34) and does not have poles at x1 = a or x1 = ū 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 Eq. (32) with f given by Eq. (35).
We obtain the remaining differentiable mapping x̄2(x1, x2) in this zone by a ”parabola-construction” as in the case of zone ④, this time, however, in the x2-directionEq. (29) except that δ is interchanged with γ and with . The ellipse parameters taken from the x̄2 transformation in zone ③ (see Eq. (26)) are
Zone ⑥: In this zone we only need to find the mapping of one component because any point in this area obeys the symmetry relationEq. (36)) to the x̄1 mapping in the area [ū, 1−ū]×[a,ū]. With the abbreviations Δ := x+ − x−, v̄:= ū−0.5 and the boundary conditions read
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 τ.
Finally, we find it very instructive to compare the transformed permittivities for the different mesh types. In Fig. 12, we display the -component of the effective permittivity, i.e., 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).
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
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.
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.
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.
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 . 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 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.
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 toEq. (44). An ellipse with center (x0, y0) that is rotated by an angle ϕ with semi-axes c and d is described by Eq. (44) with Eq. (45) yields 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. . 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. . In Fig. 17(a), we display an illustrative example of such a one-dimensional compression.
The transformation function of Ref.  is9] 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̄1 and x̄2. Furthermore, we assume that the length of the interval that is mapped onto [x̄1, x̄2] is given—it is referred to as Δx. Then, we perform the mapping Eq. (47) as if the first material interface lies at 0.Fig. 17(a)—we now have a compression at two interfaces with the correct spacing between them, albeit at the wrong locations. In order to achieve the compression in the interior in the unit cell, we first shift the function by x̄1 upwards, see Fig. 17(b). The point where the compression function intersects with the unit cell edge, denoted x̃, is Eq. (52) is a transcendental equation and, therefore, we have to find the point x̃ numerically. Next, we shift the function to the right by L̃ = L − x̃, see Fig. 17(c). The part that is still outside the unit cell (marked in red) is attached at the left side of the unit cell—this makes sense due to periodic boundaries in the FMM. Finally, in Fig. 17(d) we depict how the equidistant, vertical coordinate lines are now mapped to the compressed horizontal lines. The compression function for two material interface points x̄1 and x̄2 in the interior of the unit cell then reads 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̄2 − x̄1)/Δx which leads to a vanishing prefactor of the sine in Eq. (47). This is how we have obtained the compression in Fig. 18(b).
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 . 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.,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.
As the cylindrical fiber is centered in the square unit cell we can apply the same compression function in x̄1- and x̄2-direction. The material surface, i.e., the compression points are x̄1 = x− and x̄2 = x+ for the nondifferentiable mesh and x̄1 = α and x̄2 = 1 − α for the differentiable mesh (see section 4 and 5, respectively). In Figs. 19(b) to 19(d) we display the dependence of the error on the compression parameters G and Δx for the nondifferentiable mesh (τ = 0, panel (b)) and the differentiable mesh with parameters τ = 0.002 (panel (c)) and τ = 0.015 (panel (d)). Each computation has been performed with 997 plane waves in conjunction with a discretization in transformed space of 1024 × 1024 points. In these graphs the parameter G has been stepped in 0.01 intervals and Δx in 0.0025 intervals. We have used the same parameter stepping for G and Δx in Figs. 19 to 22. For illustrative purposes, we display in Fig. 19(d) selected compression functions for vastly different parameter values G and Δx.
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.
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. . 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 · 1016 Hz and a damping coefficient γ = 1.1536 · 1014 Hz, corresponding to gold . 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).
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.
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).
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]
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]
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. LoveOptical 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]