Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group
Journal Home About Issues in Progress Current Issue All Issues Feature Issues

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

Jens Küchenmeister, Thomas Zebrowski, and Kurt Busch

Open Access Open Access

Abstract

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

© 2012 Optical Society of America

1. Introduction

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

A significant advancement has been achieved when Lalanne and Morris [3] and Granet and Guizal [4] 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]. 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]. 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. [10] 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.

 figure: Fig. 1

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

Download Full Size | PDF

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 form

x¯1=x¯1(x1,x2),
x¯2=x¯2(x1,x2),
x¯3=x3.
Our aim is to solve Maxwell’s curl equations which are given in covariant form as [11]
ξρατ∂σEτ=ik0gμρσHσ,
ξρστ∂σHτ=−ik0gερσEσ.
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 k0 = ω/c is the vacuum wave number (c denotes the vacuum speed of light and ω the frequency). The metric tensor is
gρσ=∂xρ∂x¯τ∂xσ∂x¯τ,
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
ερσ=∂xρ∂x¯τ∂xσ∂x¯χε¯τχ,
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̄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 [2]. 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]. 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. [13].

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̄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 gερσ, with the square root of the reciprocal of the metric determinant g and ερσ like in Eq. (7). Therefore, we refer to gερσ from now on as the effective permittivity. For our present example, we depict the 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 ∂x1/∂x̄1. It can be obtained by (for other derivatives, the situation is completely analogous)

∂x1∂x¯1=(∂x¯1∂x1∂x¯2∂x2−∂x¯1∂x2∂x¯2∂x1)−1︸=|J|−1∂x¯2∂x2.
If the coordinate lines of the mesh are parallel in x̄1- and x̄2-direction at one point, the determinant of the Jacobian J of the coordinate transformation will become zero at this point. This leads to a diverging derivative ∂x1/∂x̄1 and, thereby, to a diverging effective permittivity.

 figure: Fig. 2

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̄1- and x̄2-direction is marked with a red circle. Panel (c) displays the 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.

Download Full Size | PDF

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.

 figure: Fig. 3

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 r1), a small, inner circle (red, radius r2), and two tiny circles at the crescent tips (radii r3 and r4). To define the orientation of the crescent within the unit cell, the centers (M1, N1) and (M2, N2) of the inner and outer circles need to be specified. The centers (M3, N3) and (M4, N4) 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.

Download Full Size | PDF

The crescent is composed of four circle arcs described by

(x1−Mi)2+(x2−Ni)2=ri2,             i=1,…,4,
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
N3=(r2+r3)2−(r1−r3)2+N12−N222N1−2N2,     M3=M1−(r1−r3)2−(N3−N1)2.
The point P=(xP1,xP2) where the tip circle 3 meets the outer circle 1 is given by
xP1=M3−r3cos(arctan(N3−N1M1−M3)),       xP2=N3+r32−(xP1−M3)2,
while the intersection point D=(xD1,xD2) of circle 3 and 2 is
xD1=M3+r3cos(arctan(N2−N3M2−M3)),     xD2=N3+r32−(xD1−M3)2.

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

LT(c,c¯,d,d¯,x)=d¯−c¯d−cx+c¯−cd¯−c¯d−c
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).

 figure: Fig. 4

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.

Download Full Size | PDF

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=(xR1,xR2) is given by the intersection of the (Cartesian) coordinate line xP1 with the structure surface. Then, points Q and R follow from the structure’s symmetry.

R=(xP1,N1−r12−(xR1−M1)2),      Q=(2M1−xP1,xP2),       S=(xQ1,xR2).
The (Cartesian) coordinate lines passing through any two of these characteristic points divide the unit cell into the zones ① to ⑨. Second, the so selected (Cartesian) coordinate lines have to be deformed in order to follow the crescent’s surface and this surface has to be parametrized (see Eq. (9)). Third and final, we have to find the two-dimensional transformation function (x̄1(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
x¯2(x1,x2)=x2,         (x1,x2)∈④.
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
x¯1(x1,x2)=LT(0,    0︸map left edge onto itself,xP1,      M1−r12−(x2−N1)2︸map PR coordinate line on circle arc,x1)(13)=M1−r12−(x2−N1)2x1xP1,(x1,x2)∈④.
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].

 figure: Fig. 5

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.

Download Full Size | PDF

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

CAL/R(i,x2)=Mi∓ri2−(x2−Ni)2,CAT/B(i,x1)=Ni±ri2−(x1−Mi)2,i=1,…,4.
Here, the subscripts L, R, T, and B refer to left, right, top, and bottom, respectively.

Finally, upon introducing the point E=(xE1,xE2)=(2M1−xD1,xD2), we arrive at the complete coordinate transformation as

x¯1(x1,x2)={x1(x1,x2)∈①,②,③,⑦,⑧,⑨LT(0,0,xP1,CAL(1,x2),x1)(x1,x2)∈④LT(xP1,CAL(1,x2),xQ1,CAR(1,x2),x1)(x1,x2)∈⑤LT(xQ1,CAR(1,x2),1,1,x1)(x1,x2)∈⑥.
The x̄2 mapping is more complex but the construction follows the same principles. It reads
x¯2(x1,x2)={x2(x1,x2)∈①,③​,④,⑥,⑦,⑨LT(L,L,xP2,CAT(3,x1),x2)x1≤xD1LT(L,L,xP2,CAB(2,x1),x2)x1∈(xD1,xE1)LT(L,L,xP2,CAT(4,x1),x2)xE1≤x1}(x1,x2)∈②LT(xP2,CAT(3,x1),xR2,CAB(1,x1),x2)x1≤xD1LT(xP2,CAB(2,x1),xR2,CAB(1,x2),x2)x1∈(xD1,xE1)LT(xP2,CAT(4,x1),xR2,CAB(1,x1),x2)xE1≤x1}(x1,x2)∈⑤LT(xR2,CAB(1,x1),0,0,x2)(x1,x2)∈⑧.

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.

 figure: Fig. 6

Fig. 6 Nondifferentiable meshes for a crescent with parameters r1 = 0.25, r2 = 0.2, r3 = r4 = 0.02, M1 = 0.45, N1 = 0.38, M2 = M1, N2 = 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, Δx1 = 0.4, x¯11=xP1, x¯21=xQ1, Δx2 = 0.4, x¯12=xR2, x¯22=xP2 has been applied (see section 8 for details). Panel (c) features a close-up of (a) of the left crescent tip.

Download Full Size | PDF

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.

 figure: Fig. 7

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

Download Full Size | PDF

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̄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 x±=0.5±r2.

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(x1,a,b,x−)=b(x1−a)2+x−
for the parabola g we assure continuity and differentiability at a. Demanding continuity and differentiability at ū determines the parameters a and b.
a=−2u¯−0.5(−r2−(u¯−0.5)2+0.5−x−)r2−(u¯−0.5)2+u¯,
b=14(u¯−0.5)2[(−r2−(u¯−0.5)2+0.5−x−)(r2−(u¯−0.5)2)]−1.
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 CA±(x1)=±r2−(x1−0.5)2+0.5. The mapping for the horizontal lines, i.e., the component x̄2 then reads
x¯2(x1,x2)={x2,x1∈[0,a)∪(1−a,1],x2∈[0,1]LT(0,0,x−,g(x1,a,b,x−),x2),x1∈[a,u¯]LT(0,0,x−,CA−(x1),x2)x1∈(u¯,1−u¯)LT(0,0,x−,g(x1,1−a,b,x−),x2),x1∈[1−u¯,1−a]}x2∈[0,x−]LT(x−,g(x1,a,b,x−),x+,g(x1,a,−b,x+),x2),x1∈[a,u¯]LT(x−,CA−(x1),x+,CA+(x1),x2)x1∈(u¯,1−u¯)LT(x−,g(x1,1−a,b,x−),x+,g(x1,1−a,−b,x+),x2),x1∈[1−u¯,1−a]}x2∈(x−,x+)LT(x+,g(x1,a,−b,x+),1,1,x2),x1∈[a,u¯]LT(x+,CA+(x1),1,1,x2),x1∈(u¯,1−u¯)LT(x+,g(x1,1−a,−b,x+),1,1,x2),x1∈[1−u¯,1−a]}x2∈[x+,1]

 figure: Fig. 8

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.

Download Full Size | PDF

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

 figure: Fig. 9

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.

Download Full Size | PDF

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

 figure: Fig. 10

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.

Download Full Size | PDF

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

x¯1(x1,x2)=x1,                    x¯2(x1,x2)=x2,          (x1,x2)∈①,
x¯1(x1,x2)=x1x−(0.5−r2−(x2−0.5)2),     x¯2(x1,x2)=x2,       (x1,x2)∈②
x¯1(x1,x2)=2x1−1x+−x−r2−(x2−0.5)2+12,     x¯2(x1,x2)=x¯1(x2,x1),        (x1,x2)∈③.

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

r2=(x¯1−x01)2δ+(x¯2−x02)2γ
with the parameters x01, x02, δ and γ. From this, we can infer the ellipse parameters in zone ② by comparing with the mapping in Eq. (25): x01=12x1x−, x02=0.5, δ=(x−x1)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
x¯1(x1,x2)=b(x1)(x2−a)2+x1,       x¯2(x1,x2)=x2,        (x1,x2)∈④.
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(x1) parametrizes the family of parabolas. Continuity and differentiability at x̄2 = ū require that
b(x1)=14(u¯=x02)2δγ2(r2−γ(u¯−x02)2)(−1δ(r2−γ(u¯−x02)2)+x01−x1).

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

x¯1(a,x2)(25)=ax−(0.5−r2−(x2−0.5)2),   x¯1(u¯,x2)(26)=2u¯−1x+−x−r2−(x2−0.5)2+0.5,
whereas differentiability requires
∂x¯1∂x1|(a,x2)(25)=1x−(0.5−r2−(x2−0.5)2),∂x¯1∂x1|(u¯,x2)(26)=2r2−(x2−0.5)2x+−x−,∂x¯1∂x2|(a,x2)(25)=ax−x2−0.5r2−(x2−0.5)2,∂x¯1∂x2|(u¯,x2)(26)=−2u¯−1x+−x−x2−0.5r2−(x2−0.5)2.

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

x¯1(x1,x2)=x¯1(a,x2)x1−u¯a−u¯+x¯1(u¯,x2)x1−au¯−a+f(x1,x2)(x1−u¯)(x1−a).
This ansatz fulfills the continuity requirements of Eq. (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
∂x¯1∂x1(x1,x2)=x¯1(a,x2)a−u¯+x¯1(u¯,x2)u¯−a+∂f∂x1(x1,x2)(x1−u¯)(x1−a)+f(x1,x2)(x1−a)+f(x1,x2)(x1−u¯).
Evaluating Eq. (33) at ū and a and comparing with Eq. (31) results in
f(u¯,x2)=1u¯−a(−x¯1(a,x2)a−u¯−x¯1(u¯,x2)u¯−a+∂x¯1∂x1|(u¯,x2)),f(a,x2)=1a−u¯(−x¯1(a,x2)a−u¯−x¯1(u¯,x2)u¯−a+∂x¯1∂x1|(a,x2)).
We would like to emphasize that every function f that fulfills 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
f(x1,x2)=f(u¯,x2)−f(a,x2)u¯−ax1+f(u¯,x2)−f(u¯,x2)−f(a,x2)u¯−au¯.
Thereby, we acquired a differentiable mapping x̄1(x1, x2) in zone ⑤—namely 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-direction

x¯2(x1,x2)=b(x2)(x1−a)2+x2.
Then, b(x2) is identical to Eq. (29) except that δ is interchanged with γ and x01 with x02. The ellipse parameters taken from the x̄2 transformation in zone ③ (see Eq. (26)) are
x01=0.5,     x02=0.5,     γ=(x+−x−2x2−1)2,        δ=1.

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

(x¯1(x1,x2),x¯2(x1,x2))=(x¯2(x2,x1),x¯1(x2,x1)).
The boundary conditions for this zone are derived from the mappings in zones ④ and ⑤. In particular, we use symmetry to connect the x̄2 mapping in zone ⑤ (Eq. (36)) to the x̄1 mapping in the area [ū, 1−ū]×[a,ū]. With the abbreviations Δ := x+ − x−, v̄:= ū−0.5 and s:=r2−v¯2 the boundary conditions read
x¯1(a,x2)(28)=b′(a)(x2−a)2+a,withb′(a)(29)=14ax−v¯2(−s+0.5−x−)s2,x¯1(u¯,x2)(36)=b″(u¯)(x2−a)2+u¯,withb″(u¯)(37)=12v¯3(s−0.5Δ)Δs2,
∂x¯1∂x2|(a,x2)(28)=2b′(a)(x2−a),∂x¯1∂x1|(a,x2)(28)=14v¯2(x2−a)2(−s+0.5−x−)x−s2+1,∂x¯2∂x2|(u¯,x2)(36)=2b″(u¯)(x2−a),∂x¯1∂x1|(u¯,x2)(36)=12v¯2(x2−a)2(s−0.5Δ)Δs2+1.

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

 figure: Fig. 11

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.

Download Full Size | PDF

Finally, we find it very instructive to compare the transformed permittivities for the different mesh types. In Fig. 12, we display the gε11-component of the effective permittivity, i.e., 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).

 figure: Fig. 12

Fig. 12 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.

Download Full Size | PDF

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=x¯1(α,0.5),          (α,0.5)∈⑤

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.

 figure: Fig. 13

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.

Download Full Size | PDF

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.

 figure: Fig. 14

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.

Download Full Size | PDF

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: Fig. 15

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.

Download Full Size | PDF

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

xα=xrec+yrecsinα↔xrec=xα−yαsinαcosα,yα=yreccosα↔yrec=yα1cosα.
Here, the coordinates xα, yα, and xrec, yrec 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.

 figure: Fig. 16

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.

Download Full Size | PDF

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−12(1+sinα))2+(y−12cosα)2=r2},
Circlerec(42)={(x−ysinαcosα,y1cosα):(x−12(1+sinα))2+(y−12cosα)2=r2},={(x,y):(x+ysinα−12(1+sinα))2+(ycosα−12cosα)2=r2}.
In order to perform the mesh generation analytically, we have to extract the ellipse parameters from Eq. (44). An ellipse with center (x0, y0) that is rotated by an angle ϕ with semi-axes c and d is described by
(x−x0y−y0)(cosϕ−sinϕsinϕcosϕ)(c−200d−2)(cosϕsinϕ−sinϕcosϕ)(x−x0y−y0)=1.
Comparing Eq. (44) with Eq. (45) yields
c=r(1−sinα)−1/2,   d=r(1+sinα)−1/2,      ϕ=−45°,      (x0,y0)=(0.5,0.5).
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]. 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]. In Fig. 17(a), we display an illustrative example of such a one-dimensional compression.

 figure: Fig. 17

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

Download Full Size | PDF

The transformation function of Ref. [9] is

x¯(x)=α+βx+γ2πsin(2πx−xl−1xl−xl−1),      x∈[xl−1,xl],       l=2,…,n
with
α=xlx¯l−1−xl−1x¯lxl−xl−1,      β=x¯l−x¯l−1xl−xl−1,       γ=(xl−xl−1)G−(x¯l−x¯l−1),
where {xl : l = 1,...,n} and {x̄l : l = 1,...,n} denote two sets of points. By design, this means that x̄(xl) = x̄l for l = 1,..., n. Further, G denotes the slope of the function at all positions xl. 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 [xl−1, xl] is mapped onto the interval [x̄l−1, x̄l]. In order to maintain periodicity with period L, one chooses x̄1 = 0 and x̄n = L. This means that the coordinate lines are compressed at the unit cell edges. For the structures considered in Ref. [9] 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.

x∈[0,Δx]:  x¯(x)=x¯2−x¯1Δxx+GΔx−(x¯2−x¯1)2πsin(2πxΔx),
x∈[Δx,L]:   x¯(x)=α˜+β˜x+γ˜sin(2πx−ΔxL−Δx).
Here, we have introduced the abbreviations
α˜=L(x¯2−x¯1)−ΔxLL−Δx,      β˜=L−(x¯2−x¯1)L−Δx,     γ˜=G(L−Δx)−(L−(x¯2−x¯1))2π.
The result is depicted in Fig. 17(a)—we now have a compression at two interfaces with the correct spacing between them, albeit at the wrong locations. In order to achieve the compression in the interior in the unit cell, we first shift the function by x̄1 upwards, see Fig. 17(b). The point where the compression function intersects with the unit cell edge, denoted x̃, is
L=α˜+β˜x˜+γ˜sin(2πx˜−ΔxL−Δx)+x¯1
Unfortunately, Eq. (52) is a transcendental equation and, therefore, we have to find the point x̃ numerically. Next, we shift the function to the right by L̃ = L − x̃, see Fig. 17(c). The part that is still outside the unit cell (marked in red) is attached at the left side of the unit cell—this makes sense due to periodic boundaries in the FMM. Finally, in Fig. 17(d) we depict how the equidistant, vertical coordinate lines are now mapped to the compressed horizontal lines. The compression function for two material interface points x̄1 and x̄2 in the interior of the unit cell then reads
x∈[0,L˜]:   x¯(x)=α˜+β˜(x+x˜)+γ˜sin(2πx+x˜−ΔxL−Δx)+x¯1−L,
x∈(L˜,L˜+Δx):       x¯(x)=x¯2−x¯1Δx(x−L˜)+GΔx−(x¯2−x¯2)2πsin(2πx−L˜Δx)+x¯1,
x∈[L˜+Δx,L]:   x¯(x)=α˜+β˜(x−L˜)+γ˜sin(2πx−L˜−ΔxL−Δx)+x¯1.
We want to stress that one needs to start with the separation Δx of the interfaces instead of starting with their exact locations x1 and x2 (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̄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).

 figure: Fig. 18

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̄1 and x̄2) is applied in the vertical direction. Panel (c) depicts the resulting mesh.

Download Full Size | PDF

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]. 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{|neff,analyt,1st mode−neff,num,1st mode|neff,analyt,1st mode,…,|neff,analyt,10th mode−neff,num,10th mode|neff,analyt,10th mode},
where the effective refractive index of an eigenmode with propagation constant γ is given by neff = λγ/(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.

 figure: Fig. 19

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̄1- and x̄2-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.

Download Full Size | PDF

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.

 figure: Fig. 20

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.

Download Full Size | PDF

 figure: Fig. 21

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.

Download Full Size | PDF

 figure: Fig. 22

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.

Download Full Size | PDF

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.

 figure: Fig. 23

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.

Download Full Size | PDF

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 · 1016 Hz and a damping coefficient γ = 1.1536 · 1014 Hz, corresponding to gold [15]. 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).

 figure: Fig. 24

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.

Download Full Size | PDF

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.

 figure: Fig. 25

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.

Download Full Size | PDF

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

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (25)
Fig. 1
Fig. 1 Schematic view of a typical system of interest in FMM computations. The system is finite in x̄3-direction and periodic with respect to the x̄1x̄2 plane. kin denotes the wave vector of an incident plane wave. Ox̄1x̄2x̄3 describes a Cartesian coordinate system.

View in Article | Download Full Size | PDF

Fig. 2
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̄1- and x̄2-direction is marked with a red circle. Panel (c) displays the 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.

View in Article | Download Full Size | PDF

Fig. 3
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 r1), a small, inner circle (red, radius r2), and two tiny circles at the crescent tips (radii r3 and r4). To define the orientation of the crescent within the unit cell, the centers (M1, N1) and (M2, N2) of the inner and outer circles need to be specified. The centers (M3, N3) and (M4, N4) 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.

View in Article | Download Full Size | PDF

Fig. 4
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.

View in Article | Download Full Size | PDF

Fig. 5
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.

View in Article | Download Full Size | PDF

Fig. 6
Fig. 6 Nondifferentiable meshes for a crescent with parameters r1 = 0.25, r2 = 0.2, r3 = r4 = 0.02, M1 = 0.45, N1 = 0.38, M2 = M1, N2 = 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, Δx1 = 0.4, x ¯ 1 1 = x P 1, x ¯ 2 1 = x Q 1, Δx2 = 0.4, x ¯ 1 2 = x R 2, 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.

View in Article | Download Full Size | PDF

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

View in Article | Download Full Size | PDF

Fig. 8
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.

View in Article | Download Full Size | PDF

Fig. 9
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.

View in Article | Download Full Size | PDF

Fig. 10
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.

View in Article | Download Full Size | PDF

Fig. 11
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.

View in Article | Download Full Size | PDF

Fig. 12
Fig. 12 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.

View in Article | Download Full Size | PDF

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

View in Article | Download Full Size | PDF

Fig. 14
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.

View in Article | Download Full Size | PDF

Fig. 15
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.

View in Article | Download Full Size | PDF

Fig. 16
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.

View in Article | Download Full Size | PDF

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

View in Article | Download Full Size | PDF

Fig. 18
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̄1 and x̄2) is applied in the vertical direction. Panel (c) depicts the resulting mesh.

View in Article | Download Full Size | PDF

Fig. 19
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̄1- and x̄2-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.

View in Article | Download Full Size | PDF

Fig. 20
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.

View in Article | Download Full Size | PDF

Fig. 21
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.

View in Article | Download Full Size | PDF

Fig. 22
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.

View in Article | Download Full Size | PDF

Fig. 23
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.

View in Article | Download Full Size | PDF

Fig. 24
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.

View in Article | Download Full Size | PDF

Fig. 25
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.

View in Article | Download Full Size | PDF

Equations (56)

Equations on this page are rendered with MathJax. Learn more.

x ¯ 1 = x ¯ 1 ( x 1 , x 2 ) ,
x ¯ 2 = x ¯ 2 ( x 1 , x 2 ) ,
x ¯ 3 = x 3 .
ξ ρ α τ ∂ σ E τ = i k 0 g μ ρ σ H σ ,
ξ ρ σ τ ∂ σ H τ = − i k 0 g ε ρ σ E σ .
g ρ σ = ∂ x ρ ∂ x ¯ τ ∂ x σ ∂ x ¯ τ ,
ε ρ σ = ∂ x ρ ∂ x ¯ τ ∂ x σ ∂ x ¯ χ ε ¯ τ χ ,
∂ x 1 ∂ x ¯ 1 = ( ∂ x ¯ 1 ∂ x 1 ∂ x ¯ 2 ∂ x 2 − ∂ x ¯ 1 ∂ x 2 ∂ x ¯ 2 ∂ x 1 ) − 1 ︸ = | J | − 1 ∂ x ¯ 2 ∂ x 2 .
( x 1 − M i ) 2 + ( x 2 − N i ) 2 = r i 2 ,                           i = 1 , … , 4 ,
N 3 = ( 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 − ( r 1 − r 3 ) 2 − ( N 3 − N 1 ) 2 .
x P 1 = M 3 − r 3 cos ( arctan ( N 3 − N 1 M 1 − M 3 ) ) ,               x P 2 = N 3 + r 3 2 − ( x P 1 − M 3 ) 2 ,
x D 1 = M 3 + r 3 cos ( arctan ( N 2 − N 3 M 2 − M 3 ) ) ,           x D 2 = N 3 + r 3 2 − ( x D 1 − M 3 ) 2 .
L T ( c , c ¯ , d , d ¯ , x ) = d ¯ − c ¯ d − c x + c ¯ − c d ¯ − c ¯ d − c
R = ( x P 1 , N 1 − 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 ) .
x ¯ 2 ( x 1 , x 2 ) = x 2 ,                   ( x 1 , x 2 ) ∈ ④ .
x ¯ 1 ( x 1 , x 2 ) = L T ( 0 ,         0 ︸ map   left   edge   onto   itself , x P 1 ,             M 1 − r 1 2 − ( x 2 − N 1 ) 2 ︸ map   P R   coordinate   line   on   circle   arc , x 1 ) ( 13 ) = M 1 − r 1 2 − ( x 2 − N 1 ) 2 x 1 x P 1 , ( x 1 , x 2 ) ∈ ④ .
C A L / R ( i , x 2 ) = M i ∓ r i 2 − ( x 2 − N i ) 2 , C A T / B ( i , x 1 ) = N i ± r i 2 − ( x 1 − M i ) 2 , i = 1 , … , 4.
x ¯ 1 ( x 1 , x 2 ) = { x 1 ( x 1 , x 2 ) ∈ ① , ② , ③ , ⑦ , ⑧ , ⑨ L T ( 0 , 0 , x P 1 , C A L ( 1 , x 2 ) , x 1 ) ( x 1 , x 2 ) ∈ ④ 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 ) ∈ ⑤ L T ( x Q 1 , C A R ( 1 , x 2 ) , 1 , 1 , x 1 ) ( x 1 , x 2 ) ∈ ⑥ .
x ¯ 2 ( x 1 , x 2 ) = { x 2 ( x 1 , x 2 ) ∈ ① , ③ ​ , ④ , ⑥ , ⑦ , ⑨ L T ( L , L , x P 2 , C A T ( 3 , x 1 ) , x 2 ) x 1 ≤ x D 1 L T ( L , L , x P 2 , C A B ( 2 , x 1 ) , x 2 ) x 1 ∈ ( x D 1 , x E 1 ) L T ( L , L , x P 2 , C A T ( 4 , x 1 ) , x 2 ) x E 1 ≤ x 1 } ( x 1 , x 2 ) ∈ ② L T ( x P 2 , C A T ( 3 , x 1 ) , x R 2 , C A B ( 1 , x 1 ) , x 2 ) x 1 ≤ 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 ∈ ( 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 ≤ x 1 } ( x 1 , x 2 ) ∈ ⑤ L T ( x R 2 , C A B ( 1 , x 1 ) , 0 , 0 , x 2 ) ( x 1 , x 2 ) ∈ ⑧ .
g ( x 1 , a , b , x − ) = b ( x 1 − a ) 2 + x −
a = − 2 u ¯ − 0.5 ( − r 2 − ( u ¯ − 0.5 ) 2 + 0.5 − x − ) r 2 − ( u ¯ − 0.5 ) 2 + u ¯ ,
b = 1 4 ( u ¯ − 0.5 ) 2 [ ( − r 2 − ( u ¯ − 0.5 ) 2 + 0.5 − x − ) ( r 2 − ( u ¯ − 0.5 ) 2 ) ] − 1 .
x ¯ 2 ( x 1 , x 2 ) = { x 2 , x 1 ∈ [ 0 , a ) ∪ ( 1 − a , 1 ] , x 2 ∈ [ 0 , 1 ] L T ( 0 , 0 , x − , g ( x 1 , a , b , x − ) , x 2 ) , x 1 ∈ [ a , u ¯ ] L T ( 0 , 0 , x − , C A − ( x 1 ) , x 2 ) x 1 ∈ ( u ¯ , 1 − u ¯ ) L T ( 0 , 0 , x − , g ( x 1 , 1 − a , b , x − ) , x 2 ) , x 1 ∈ [ 1 − u ¯ , 1 − a ] } x 2 ∈ [ 0 , x − ] L T ( x − , g ( x 1 , a , b , x − ) , x + , g ( x 1 , a , − b , x + ) , x 2 ) , x 1 ∈ [ a , u ¯ ] L T ( x − , C A − ( x 1 ) , x + , C A + ( x 1 ) , x 2 ) x 1 ∈ ( u ¯ , 1 − u ¯ ) L T ( x − , g ( x 1 , 1 − a , b , x − ) , x + , g ( x 1 , 1 − a , − b , x + ) , x 2 ) , x 1 ∈ [ 1 − u ¯ , 1 − a ] } x 2 ∈ ( x − , x + ) L T ( x + , g ( x 1 , a , − b , x + ) , 1 , 1 , x 2 ) , x 1 ∈ [ a , u ¯ ] L T ( x + , C A + ( x 1 ) , 1 , 1 , x 2 ) , x 1 ∈ ( u ¯ , 1 − u ¯ ) L T ( x + , g ( x 1 , 1 − a , − b , x + ) , 1 , 1 , x 2 ) , x 1 ∈ [ 1 − u ¯ , 1 − a ] } x 2 ∈ [ x + , 1 ]
x ¯ 1 ( x 1 , x 2 ) = x 1 ,                                         x ¯ 2 ( x 1 , x 2 ) = x 2 ,                     ( x 1 , x 2 ) ∈ ① ,
x ¯ 1 ( x 1 , x 2 ) = x 1 x − ( 0.5 − r 2 − ( x 2 − 0.5 ) 2 ) ,           x ¯ 2 ( x 1 , x 2 ) = x 2 ,               ( x 1 , x 2 ) ∈ ②
x ¯ 1 ( x 1 , x 2 ) = 2 x 1 − 1 x + − x − r 2 − ( x 2 − 0.5 ) 2 + 1 2 ,           x ¯ 2 ( x 1 , x 2 ) = x ¯ 1 ( x 2 , x 1 ) ,                 ( x 1 , x 2 ) ∈ ③ .
r 2 = ( x ¯ 1 − x 0 1 ) 2 δ + ( x ¯ 2 − x 0 2 ) 2 γ
x ¯ 1 ( x 1 , x 2 ) = b ( x 1 ) ( x 2 − a ) 2 + x 1 ,               x ¯ 2 ( x 1 , x 2 ) = x 2 ,                 ( x 1 , x 2 ) ∈ ④ .
b ( x 1 ) = 1 4 ( u ¯ = x 0 2 ) 2 δ γ 2 ( r 2 − γ ( u ¯ − x 0 2 ) 2 ) ( − 1 δ ( r 2 − γ ( u ¯ − x 0 2 ) 2 ) + x 0 1 − x 1 ) .
x ¯ 1 ( a , x 2 ) ( 25 ) = a x − ( 0.5 − r 2 − ( x 2 − 0.5 ) 2 ) ,       x ¯ 1 ( u ¯ , x 2 ) ( 26 ) = 2 u ¯ − 1 x + − x − r 2 − ( x 2 − 0.5 ) 2 + 0.5 ,
∂ x ¯ 1 ∂ x 1 | ( a , x 2 ) ( 25 ) = 1 x − ( 0.5 − r 2 − ( x 2 − 0.5 ) 2 ) , ∂ x ¯ 1 ∂ x 1 | ( u ¯ , x 2 ) ( 26 ) = 2 r 2 − ( x 2 − 0.5 ) 2 x + − x − , ∂ x ¯ 1 ∂ x 2 | ( a , x 2 ) ( 25 ) = a x − x 2 − 0.5 r 2 − ( x 2 − 0.5 ) 2 , ∂ x ¯ 1 ∂ x 2 | ( u ¯ , x 2 ) ( 26 ) = − 2 u ¯ − 1 x + − x − x 2 − 0.5 r 2 − ( x 2 − 0.5 ) 2 .
x ¯ 1 ( x 1 , x 2 ) = x ¯ 1 ( a , x 2 ) x 1 − u ¯ a − u ¯ + x ¯ 1 ( u ¯ , x 2 ) x 1 − a u ¯ − a + f ( x 1 , x 2 ) ( x 1 − u ¯ ) ( x 1 − a ) .
∂ x ¯ 1 ∂ x 1 ( x 1 , x 2 ) = x ¯ 1 ( a , x 2 ) a − u ¯ + x ¯ 1 ( u ¯ , x 2 ) u ¯ − a + ∂ f ∂ x 1 ( x 1 , x 2 ) ( x 1 − u ¯ ) ( x 1 − a ) + f ( x 1 , x 2 ) ( x 1 − a ) + f ( x 1 , x 2 ) ( x 1 − u ¯ ) .
f ( u ¯ , x 2 ) = 1 u ¯ − a ( − x ¯ 1 ( a , x 2 ) a − u ¯ − x ¯ 1 ( u ¯ , x 2 ) u ¯ − a + ∂ x ¯ 1 ∂ x 1 | ( u ¯ , x 2 ) ) , f ( a , x 2 ) = 1 a − u ¯ ( − x ¯ 1 ( a , x 2 ) a − u ¯ − x ¯ 1 ( u ¯ , x 2 ) u ¯ − a + ∂ x ¯ 1 ∂ x 1 | ( a , x 2 ) ) .
f ( x 1 , x 2 ) = f ( u ¯ , x 2 ) − f ( a , x 2 ) u ¯ − a x 1 + f ( u ¯ , x 2 ) − f ( u ¯ , x 2 ) − f ( a , x 2 ) u ¯ − a u ¯ .
x ¯ 2 ( x 1 , x 2 ) = b ( x 2 ) ( x 1 − a ) 2 + x 2 .
x 0 1 = 0.5 ,           x 0 2 = 0.5 ,           γ = ( x + − x − 2 x 2 − 1 ) 2 ,                 δ = 1.
( x ¯ 1 ( x 1 , x 2 ) , x ¯ 2 ( x 1 , x 2 ) ) = ( x ¯ 2 ( x 2 , x 1 ) , x ¯ 1 ( x 2 , x 1 ) ) .
x ¯ 1 ( a , x 2 ) ( 28 ) = b ′ ( a ) ( x 2 − a ) 2 + a , with b ′ ( a ) ( 29 ) = 1 4 a x − v ¯ 2 ( − s + 0.5 − x − ) s 2 , x ¯ 1 ( u ¯ , x 2 ) ( 36 ) = b ″ ( u ¯ ) ( x 2 − a ) 2 + u ¯ , with b ″ ( u ¯ ) ( 37 ) = 1 2 v ¯ 3 ( s − 0.5 Δ ) Δ s 2 ,
∂ x ¯ 1 ∂ x 2 | ( a , x 2 ) ( 28 ) = 2 b ′ ( a ) ( x 2 − a ) , ∂ x ¯ 1 ∂ x 1 | ( a , x 2 ) ( 28 ) = 1 4 v ¯ 2 ( x 2 − a ) 2 ( − s + 0.5 − x − ) x − s 2 + 1 , ∂ x ¯ 2 ∂ x 2 | ( u ¯ , x 2 ) ( 36 ) = 2 b ″ ( u ¯ ) ( x 2 − a ) , ∂ x ¯ 1 ∂ x 1 | ( u ¯ , x 2 ) ( 36 ) = 1 2 v ¯ 2 ( x 2 − a ) 2 ( s − 0.5 Δ ) Δ s 2 + 1.
0.5 − r = x ¯ 1 ( α , 0.5 ) ,                     ( α , 0.5 ) ∈ ⑤
x α = x rec + y rec sin α ↔ x rec = x α − y α sin α cos α , y α = y rec cos α ↔ y rec = y α 1 cos α .
Circle α = { ( x , y ) : ( x − 1 2 ( 1 + sin α ) ) 2 + ( y − 1 2 cos α ) 2 = r 2 } ,
Circle rec ( 42 ) = { ( x − y sin α cos α , y 1 cos α ) : ( x − 1 2 ( 1 + sin α ) ) 2 + ( y − 1 2 cos α ) 2 = r 2 } , = { ( x , y ) : ( x + y sin α − 1 2 ( 1 + sin α ) ) 2 + ( y cos α − 1 2 cos α ) 2 = r 2 } .
( x − x 0 y − y 0 ) ( cos ϕ − sin ϕ sin ϕ cos ϕ ) ( c − 2 0 0 d − 2 ) ( cos ϕ sin ϕ − sin ϕ cos ϕ ) ( x − x 0 y − y 0 ) = 1.
c = r ( 1 − sin α ) − 1 / 2 ,       d = r ( 1 + sin α ) − 1 / 2 ,             ϕ = − 45 ° ,             ( x 0 , y 0 ) = ( 0.5 , 0.5 ) .
x ¯ ( x ) = α + β x + γ 2 π sin ( 2 π x − x l − 1 x l − x l − 1 ) ,             x ∈ [ x l − 1 , x l ] ,               l = 2 , … , n
α = x l x ¯ l − 1 − x l − 1 x ¯ l x l − x l − 1 ,             β = x ¯ l − x ¯ l − 1 x l − x l − 1 ,               γ = ( x l − x l − 1 ) G − ( x ¯ l − x ¯ l − 1 ) ,
x ∈ [ 0 , Δ x ] :     x ¯ ( x ) = x ¯ 2 − x ¯ 1 Δ x x + G Δ x − ( x ¯ 2 − x ¯ 1 ) 2 π sin ( 2 π x Δ x ) ,
x ∈ [ Δ x , L ] :       x ¯ ( x ) = α ˜ + β ˜ x + γ ˜ sin ( 2 π x − Δ x L − Δ x ) .
α ˜ = L ( x ¯ 2 − x ¯ 1 ) − Δ x L L − Δ x ,             β ˜ = L − ( x ¯ 2 − x ¯ 1 ) L − Δ x ,           γ ˜ = G ( L − Δ x ) − ( L − ( x ¯ 2 − x ¯ 1 ) ) 2 π .
L = α ˜ + β ˜ x ˜ + γ ˜ sin ( 2 π x ˜ − Δ x L − Δ x ) + x ¯ 1
x ∈ [ 0 , L ˜ ] :       x ¯ ( x ) = α ˜ + β ˜ ( x + x ˜ ) + γ ˜ sin ( 2 π x + x ˜ − Δ x L − Δ x ) + x ¯ 1 − L ,
x ∈ ( L ˜ , L ˜ + Δ x ) :               x ¯ ( x ) = x ¯ 2 − x ¯ 1 Δ x ( x − L ˜ ) + G Δ x − ( x ¯ 2 − x ¯ 2 ) 2 π sin ( 2 π x − L ˜ Δ x ) + x ¯ 1 ,
x ∈ [ L ˜ + Δ x , L ] :       x ¯ ( x ) = α ˜ + β ˜ ( x − L ˜ ) + γ ˜ sin ( 2 π x − L ˜ − Δ x L − Δ x ) + x ¯ 1 .
error = max { | n eff , analyt , 1 st   mode − n eff , num , 1 st   mode | n eff , analyt , 1 st   mode , … , | n eff , analyt , 10 th   mode − n eff , num , 10 th   mode | n eff , analyt , 10 th   mode } ,
Select as filters
Select Topics Cancel
Top
       
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.
Privacy | Terms of Use