Matched coordinates and adaptive spatial resolution in the Fourier modal method

Thomas Weiss; Gérard Granet; Nikolay A. Gippius; Sergei G. Tikhodeev; Harald Giessen

doi:10.1364/OE.17.008051

1. Introduction

The calculation of two-dimensional gratings occurs in many parts of modern optics, such as photonic crystals, diffractive optics, and metamaterials. A natural choice for periodic problems is the Fourier modal method (FMM) or rigorous coupled wave analysis (RCWA). It has the advantage that the structure of interest can be calculated without discretization of derivatives and no approximation by meshes is needed in at least one dimension.

The main idea of the classical Fourier modal method is to divide a three-dimensional structure into several planar periodic slabs that are invariant in their normal direction (see Fig. 1(a)) and solve Maxwell’s equations in each slab by decomposing into Floquet-Bloch modes. Then, the slabs are combined under consideration of the boundary conditions. Finally, a scattering matrix formalism [1,2] is used to derive the optical properties of the whole structure.

Following from a straight-forward calculus, the only source for inaccuracies is the truncation of the involved Floquet-Fourier expansion for numerical calculations. However, if the calculated nanostructured material exhibits large changes in the permittivity or permeability function, these inaccuracies evolve into poor convergence.

As shown by Li [3], one reason for the slow convergence is that the convolutions of two functions with concurrent jump discontinuities must lead to wrong results. He formulated rules to avoid these operations and to improve the convergence. Still, the oscillating behavior of the Floquet-Fourier expansion of discontinuous functions, also known as the Gibbs phenomenon, hampers a fast convergence, especially for metal-dielectric structures. Adaptive spatial resolution (ASR) has been suggested [4,5] so that the resolution around the jump discontinuities is increased. Consequently, not merely the discontinuous function, but rather a product with partial derivatives of the coordinate transformation has to be Fourier-transformed. This product exhibits a smaller jump height and, therefore, converges faster.

Several approaches have been presented to apply the factorization rules correctly for arbitrary 2D structures [6, 7, 8, 9]. We present here an approach to include also adaptive spatial resolution using a covariant description and matched coordinates aligned along the interfaces of the structure and illustrate the method in the example of a cylinder. Furthermore, we discuss a symmetric formulation to Fourier-factorize the permittivity and permeability tensors in Maxwell’s equations correctly.

2. Coordinate transformations

A crucial part of this work is the application of coordinate transformations in Maxwell’s equations. Therefore, we describe briefly coordinate transformations using a covariant notation.

By choosing a center point, any other point in the three-dimensional space can be displayed as a coordinate vector r that is determined by the distance and the direction with respect to the center. In a Cartesian system,

r (x^{1}, x^{2}, x^{3}) = x^{m} e_{m}, m = 1,2,3 .

Note that we use the sum convention, in which a sum is carried out when an upper index and a lower index are displayed by the same letters. Owing to the orthonormality of the Cartesian coordinate system, the coordinates x^m are the projection of r to the vectors e _m. Alternatively, r can be displayed in dependence of new coordinates ̅xm:

r ({\overset{̅}{x}}^{1}, {\overset{̅}{x}}^{2}, {\overset{̅}{x}}^{3}) = x^{m} ({\overset{̅}{x}}^{1}, {\overset{̅}{x}}^{2}, {\overset{̅}{x}}^{3}) e_{m} .

If the transformation is continuous, invertible, and has no singular points, then the above equation defines the coordinate curves r(x̄¹, const, const), r(const,x̄², const), and r(const, const, x̄³) with the tangential vectors em, as well as the coordinate surfaces x̄^m(x ¹,x ²,x ³) = const with the normal vectors e̅^m. These vectors are defined as

{\overset{̅}{e}}_{m} = \frac{\partial r}{\partial {\overset{̅}{x}}^{m}}, {\overset{̅}{e}}^{m} = \nabla {\overset{̅}{x}}^{m} .

The orientation of these vectors and curves for a two-dimensional transformation is displayed in Fig. 1(b).

Fig. 1. (a) Schematic drawing of a standard structure and definition of a Cartesian coordinate system including the incidence parameters. (b) Sketch of lines with constant x̄¹ and x̄² in a plane normal to x̄³ = x ³ with the directions of the reciprocal bases.

Download Full Size | PDF

It is possible to show that the vectors in Eq. (3) form a set reciprocal bases [10] with ē^m · ē_n = δ_n^m. The representation of any vector A in one of these bases is

A = A^{m} {\overset{̅}{e}}_{m} = A_{m} {\overset{̅}{e}}^{m} .

The covariant component A_m is the projection to ē^m, while the contravariant component A^m is the projection to ē^m. This can be deduced from Eq. (4) by scalar multiplication with either ē_n or ēⁿ. In addition, we can derive:

A_{m} = ({\overset{̅}{e}}_{m} \cdot {\overset{̅}{e}}_{n}) A^{n} = \sum_{p} \frac{\partial x^{p}}{\partial {\overset{̅}{x}}^{m}} \frac{\partial x^{p}}{\partial {\overset{̅}{x}}^{n}} A_{n} = g_{mn} A^{n},

A^{m} = ({\overset{̅}{e}}^{m} \cdot {\overset{̅}{e}}^{n}) A_{n} = \sum_{p} \frac{\partial {\overset{̅}{x}}^{m}}{\partial x^{p}} \frac{\partial {\overset{̅}{x}}^{n}}{\partial x^{p}} A_{n} = g^{mn} A_{n} .

The quantities g_mn and g_mn are called metric components and conjugate metric components, respectively. They obey the relation g_mn g^np = δ_m^p.

3. Covariant formulation of Maxwell’s equations

The curl Maxwell equations can be formulated invariant under coordinate transformations [11]. In CGS units,

ξ^{mnp} \partial_{n} E_{p} = i k_{0} {\overset{̅}{μ}}^{mn} H_{n},

ξ^{mnp} \partial_{n} E_{p} = - i k_{0} {\overset{̅}{ε}}^{mn} E_{n},

with E_p and H_p being the p-covariant component of the electric or magnetic field, respectively, k ₀ = ω/c, ξ^mnp is the Levi-Civita symbol, and the ansatz exp(-iωt) has been used for the time dependence. This formulation is valid in every well-behaved coordinate system if the permittivity and permeability tensors are defined as

{\overset{̅}{ε}}^{mn} = \sqrt{g} Λ_{p}^{m} Λ_{q}^{n} ε^{pq},

{\overset{̅}{μ}}^{mn} = \sqrt{g} Λ_{p}^{m} Λ_{q}^{n} μ^{pq},

with the determinant g of the metric components, the Cartesian permittivity and permeability elements ε^pq and μ^pq, and the transformation matrices

Λ_{n}^{m} = \frac{{\partial \overset{̅}{x}}^{m}}{\partial x^{n}} .

In the special case that the permittivity or permeability is isotropic, then, for instance, ε^pq=εδ^pq, and ε̄^mn = √g g^mnε. Note that both the permittivity and permeability can depend on the coordinates x^m.

Because Eq. (7) and Eq. (8) are very general relations, we will replace ε̄^mn and μ̄^mn by ε ^mn and μ ^mn and keep in mind that the tensorial form can either be due to anisotropic media or the consequence of a coordinate transformation or both.

Let us consider planar slabs that are invariant in the normal direction x ³ of the slabs and periodic in the directions x ¹ and x ² with periods P ₁ and P ₂. In the following the lines of constant coordinates should be matched to the interfaces inside the unit cell and obey the same translation symmetry as the original problem without coordinate transformations. For the sake of simplicity we consider only rectangular lattice configurations in this work, but this is not a limitation of the method of matched coordinates.

If the permittivity and permeability tensors are of the form

\underset{=}{ε} = (\begin{array}{c} ε^{11} & ε^{12} & 0 \\ ε^{21} & ε^{22} & 0 \\ 0 & 0 & ε^{33} \end{array}), \underset{=}{μ} = (\begin{array}{c} μ^{11} & μ^{12} & 0 \\ μ^{21} & μ^{22} & 0 \\ 0 & 0 & μ^{33} \end{array}),

then we obtain by eliminating the third components of the fields in Eq. (7) and Eq. (8):

- k_{0}^{2} \partial_{3}^{2} (\begin{array}{c} E_{1} \\ E_{2} \end{array}) = 𝓛_{EH} 𝓛_{HE} (\begin{array}{c} E_{1} \\ E_{2} \end{array}),

with

𝓛_{EH} = (\begin{array}{c} - \partial_{1} {(ε^{33})}^{- 1} \partial_{2} + k_{0}^{2} μ^{21} & \partial_{1} {(ε^{33})}^{- 1} \partial_{1} + k_{0}^{2} μ^{22} \\ - \partial_{2} {(ε^{33})}^{- 1} \partial_{2} + k_{0}^{2} μ^{11} & \partial_{2} {(ε^{33})}^{- 1} \partial_{1} + k_{0}^{2} μ^{12} \end{array}),

𝓛_{HE} = (\begin{array}{c} \partial_{1} {(μ^{33})}^{- 1} \partial_{2} + k_{0}^{2} ε^{21} & - \partial_{1} {(μ^{33})}^{- 1} \partial_{1} + k_{0}^{2} ε^{22} \\ \partial_{2} {(μ^{33})}^{- 1} \partial_{2} + k_{0}^{2} ε^{11} & - \partial_{2} {(μ^{33})}^{- 1} \partial_{1} + k_{0}^{2} ε^{12} \end{array}) .

Using a decomposition in the basis of Floquet-Bloch waves exp(ik_mx^m) with the in-plane momentum defined by the incidence parameters (see Fig. 1(a)) and the periods as

k_{1} = k \sin (θ) \cos (ϕ) + \frac{2 π}{P_{1}} j_{1}, k_{2} = k \sin (θ) \sin (ϕ) + \frac{2 π}{P_{2}} j_{2}, j_{1,2} ∊ ℤ,

with k as the length of the incidence k-vector, Eq. (13) becomes a linear eigenvalue problem for k ₃ ². By limiting the total amount of Floquet-Bloch waves in a numerical calculation to a truncation order N_G, the standard tools of linear algebra can be used to solve the above problem.

4. Factorization rules

Multiplications in normal space are replaced in Fourier space by convolutions. The Floquet-Fourier coefficients of a given function f = g h are then

f_{m} = \sum_{n} g_{m - n} h_{n} .

In numerical calculations this expansion will be truncated. In comparison with the direct Floquet-Fourier expansion of f the inaccuracies caused by the truncation are increased, because the expansion of both g and h can converge very slowly. Furthermore, Li showed [3] that in the case that g and h have concurrent jump discontinuities the convolution does not converge at all. However, if f is continuous, this behavior can be avoided by considering the Floquet-Fourier expansion of h = 1/gf. The inversion of this convolution leads to the so-called inverse rule for the convolution of f = gh. Let us define the amount of coefficients f_n as a vector [f] and g _n-m as a matrix 〚g〛. With these definitions, the inverse rule is

[f] = {〚 \frac{1}{g} 〛}^{- 1} [h] .

Note that in this formulation the convergence is still hampered by the discontinuities of g and h, but the convolution converges to the correct result.

The inverse rule can be applied also on a sum f = g ₁ h ₁ + g ₂ h ₂ if f and, for instance, h ₂ are continuous. Then,

[f] = {〚 \frac{1}{g} 〛}^{- 1} [h_{1}] + {〚 \frac{1}{g_{1}} 〛}^{- 1} 〚 \frac{1}{g_{1}} g_{2} 〛 [h_{2}] .

5. Application of factorization rules to Maxwell’s equations

Factorization rules are of great importance for considering the boundary conditions at interfaces in Maxwell’s equations correctly. In a covariant formulation,Dⁿ, E^m, Bⁿ, and H^m are continuous in direction xⁿ for m ≠ n. If the permeability and permittivity are diagonal matrices, then the procedures described in [5, 12] can be applied for two-dimensional structures.

The generalized factorization rules for nonorthogonal coordinate systems and anisotropic media discussed by Li [13] lead in their straight-forward formulation to asymmetric results for totally symmetric structures. The author states that this spurious behavior vanishes for higher truncation orders or averaging of matrix elements of his equations (20) and (39). In the case of the latter the energy conservation is violated for small truncation orders. Equations (20) and (39) of [13] have been obtained by applying the correct Fourier factorization to the permittivity and permeability tensors first in one direction and then in the other direction. Other formulations can be obtained by taking some elements of Eq. (20) and some of Eq. (39). We will now pick up this idea to improve the symmetry properties and discuss the derivation briefly.

Without the loss of generality, we focus on the displacement. Because the third component of the electric field is continuous in x ¹ and x ², the third component of the displacement can be factorized by the normal rule:

{[D^{3}]}_{12} = {〚 ε^{33} 〛}_{12} {[E_{3}]}_{12} .

The indices behind the brackets indicate that the Fourier expansion has been carried out in direction x ¹ and x ², respectively.

The other components obey the following equations:

D^{1} = ε^{11} E_{1} + ε^{12} E_{2},

D^{2} = ε^{21} E_{1} + ε^{22} E_{2},

E_{1} = Δ_{ε}^{- 1} (ε^{22} D^{1} - ε^{12} D^{2}),

E_{2} = Δ_{ε}^{- 1} (- ε^{21} D^{1} + ε^{11} D^{2}),

with Δ_ε = ε ¹¹ ε ²² -ε ¹² ε ²¹.

First, let us consider Eq. (22). Here, D ² as the projection to the normal vector and E ₁ as the projection to the tangential vector (see Fig. 1(b)) are continuous with respect to x ². Therefore, we must use Eq. (19) to expand Eq. (22) in this direction:

{[D^{2}]}_{2} = {〚 \frac{1}{ε^{22}} 〛}_{2}^{- 1} {〚 \frac{1}{ε^{22}} ε^{21} 〛}_{2} {[E_{1}]}_{2} + {〚 \frac{1}{ε^{22}} 〛}_{2}^{- 1} {[E_{2}]}_{2} .

Let us now Fourier-transform Eq. (23) with respect to x ² using that E ₁ and D ² are continuous in this direction. Then,

{[D^{1}]}_{2} = {〚 \frac{1}{ε^{22}} Δ_{ε} 〛}_{2} {[E_{1}]}_{2} + {[\frac{1}{ε^{22}} ε^{12}]}_{2} {[D^{2}]}_{2 .}

By employing Eq. (25) in Eq. (26) and by using a similar procedure for Eq. (21) and Eq. (24) we derive:

{[D^{1}]}_{2} = {〈 ε^{11} 〉}_{2} {[E_{1}]}_{2} + {〈 ε^{12} 〉}_{2} {[E_{2}]}_{2},

{[D^{2}]}_{1} = {〈 ε^{21} 〉}_{1} {[E_{1}]}_{1} + {〈 ε^{22} 〉}_{1} {[E_{2}]}_{1},

with

{〈 ε^{11} 〉}_{2} = {〚 \frac{1}{ε^{22}} Δ_{ε} 〛}_{2} + {〚 \frac{1}{ε^{22}} ε^{12} 〛}_{2} {〚 \frac{1}{ε^{22}} 〛}_{2}^{- 1} {〚 \frac{1}{ε^{22}} ε^{21} 〛}_{2},

{〈 ε^{12} 〉}_{2} = {〚 \frac{1}{ε^{22}} ε^{12} 〛}_{2} {〚 \frac{1}{ε^{22}} 〛}_{2}^{- 1},

{〈 ε^{21} 〉}_{1} = {〚 \frac{1}{ε^{11}} ε^{21} 〛}_{1} {〚 \frac{1}{ε^{11}} 〛}_{1}^{- 1},

{〈 ε^{22} 〉}_{1} = {〚 \frac{1}{ε^{11}} Δ_{ε} 〛}_{1} + {〚 \frac{1}{ε^{11}} ε^{21} 〛}_{1} {〚 \frac{1}{ε^{11}} 〛}_{1}^{- 1} {〚 \frac{1}{ε^{11}} ε^{12} 〛}_{1} .

The notation 〈ε^mn〉_p denotes that the corresponding term evolves from the correct Fourier factorization of ε^mn in the direction x^p. The application of Eq. (19) for the direction x ₁ in Eq. (27) and x ² in Eq. (28) completes the application of the generalized factorization rules to the displacement:

{〈 ε^{11} 〉}_{21} = {〚 {〈 ε^{11} 〉}_{2}^{- 1} 〛}_{1}^{- 1},

{〈 ε^{12} 〉}_{21} = {〚 {〈 ε^{11} 〉}_{2}^{- 1} 〛}_{1}^{- 1} {〚 {〈 ε^{11} 〉}_{2}^{- 1} {〈 ε^{12} 〉}_{2} 〛}_{1},

{〈 ε^{21} 〉}_{12} = {〚 {〈 ε^{22} 〉}_{1}^{- 1} 〛}_{2}^{- 1} {〚 {〈 ε^{22} 〉}_{1}^{- 1} {〈 ε^{21} 〉}_{1} 〛}_{2},

{〈 ε^{22} 〉}_{12} = {〚 {〈 ε^{22} 〉}_{1}^{- 1} 〛}_{2}^{- 1} .

A similar relation can be obtained for the magnetic induction B, the permittivity μ, and the magnetic field H. Note that Eq. (33) and Eq. (36) lead to 〚〚ε ²²〛^-1 ₁〚^-1 ₂ and 〚〚ε ¹¹〛^-1 ₂〚^-1 ₁ for ε ¹² = ε ²¹ = 0, while the standard derivation (20) given in [13] results in 〚〚ε ²²〛^-1 ₁〚^-1 ₂ and 〚〚1/ε ¹¹〛^-1 ₁〚2. If ε ¹¹ = ε ²² and the unit cell is symmetric, the spectra under normal incidence depend for the formulas in [13] on the incidence polarization, although they should not. The reason is that inverting and Fourier-transforming are not commutative in a truncated Fourier space. Such behavior is avoided by using the alternative relations (33) to (36). It should be mentioned that another symmetric form to apply factorization rules can be found that converges to 〚〚1/ε ²²〛^-1 ₂〛1 and 〚〚1/ε ¹¹〛^-1 ₁〛2 for ε ¹² = ε ²¹ = 0, but the result is very complex and contains a large amount of matrix multiplications.

6. Matched coordinates and adaptive spatial resolution

The above relations describe the correct application of Fourier factorization if the coordinates match the structure so that the interfaces are aligned along the coordinate lines. Therefore, an appropriate coordinate transformation has to be found. Owing to a mismatch of the geometry of the unit cell and an arbitrary shape inside, this problem is non-trivial, especially for finding transformations with smooth partial derivatives that converge usually faster than discontinuous functions. On the other hand the matched coordinate system does not require any coordinate lines to be normal at the interfaces.

In the following we explain a possible transformation for an array of cylinders with the radius R in a square lattice configuration with period P (see Fig. 2(a)). This example can be adapted easily to other systems. Let x^m be the Cartesian coordinates and x͂^m be the matched coordinates. We can choose the coordinate system in such a way that the cylinders are centered in the unit cells and consider one unit cell defined by the lines at x^m = 0 and x^m = P. Then, the four curves

x_{\pm}^{m} ({\tilde{x}}^{n \neq m}, {\tilde{x}}_{\pm}^{m}) = {\begin{array}{c} \frac{P}{2} \pm \sqrt{R^{2} - {({\tilde{x}}^{n} - \frac{P}{2})}^{2}} for {\tilde{x}}^{n} \leq {\tilde{x}}^{n} \leq {\tilde{x}}_{+}^{n} \\ {\tilde{x}}_{\pm}^{m} elsewhere \end{array}

with m,n= 1,2 contain the interface and define lines of constant coordinates with the interface positions in the matched coordinates at x͂^m,n _± = (P ± √2R)/2. All coordinate lines in-between can be obtained by a linear change of the profile:

x^{m} ({\tilde{x}}^{n \neq m}, {\tilde{x}}^{m}) = {\begin{array}{c} \frac{{\tilde{x}}^{m}}{{\tilde{x}}_{-}^{m}} x_{-}^{m} for 0 \leq {\tilde{x}}^{m} < {\tilde{x}}_{-}^{m} \\ \frac{{\tilde{x}}^{m} - {\tilde{x}}_{-}^{m}}{{\tilde{x}}_{+}^{m} - {\tilde{x}}_{-}^{m}} x_{+}^{m} ({\tilde{x}}^{n}) + \frac{{\tilde{x}}^{m} - {\tilde{x}}_{+}^{m}}{{\tilde{x}}_{-}^{m} - {\tilde{x}}_{+}^{m}} x_{-}^{m} ({\tilde{x}}^{n}) \\ \frac{{\tilde{x}}^{m} - P}{{\tilde{x}}_{+}^{m} - P} x_{+}^{m} ({\tilde{x}}^{n}) + \frac{{\tilde{x}}^{m} - {\tilde{x}}_{+}^{m}}{P - {\tilde{x}}_{+}^{m}} P for {\tilde{x}}_{+}^{m} < {\tilde{x}}^{m} \leq P \end{array} for {\tilde{x}}_{-}^{m} \leq {\tilde{x}}^{m} \leq {\tilde{x}}_{+}^{m} .

This transformation has discontinuous partial derivatives. However, it has been shown that adaptive spatial resolution reduces the influence of discontinuities in the case of orthogonal coordinate lines [5]. The scheme can be easily combined with the matched coordinate system.

Let us use the invariance of Maxwell’s equations (7) and (8) in the covariant formalism and consider the permittivity and permeability tensors in Eq. (12) obtained in the matched coordinates. In the corresponding coordinate system with the coordinates x͂^m all interfaces are aligned along coordinate surfaces with x͂^m = const. Hence, the same scheme as described in [5] can be applied to include adaptive spatial resolution using new coordinates x̄¹ = x̄¹(x͂¹), x̄² = x̄²(x͂²), and x̄³ = x͂³ with an increased resolution at the interfaces. Then, the new permittivity tensor is

\underset{=}{\overset{̅}{ε}} = (\begin{array}{c} \frac{\partial {\tilde{x}}^{2}}{\partial {\overset{̅}{x}}^{2}} {(\frac{\partial {\tilde{x}}^{1}}{\partial {\overset{̅}{x}}^{1}})}^{- 1} ε^{11} & ε^{12} & 0 \\ ε^{21} & \frac{\partial {\tilde{x}}^{1}}{\partial {\overset{̅}{x}}^{1}} {(\frac{\partial {\tilde{x}}^{2}}{\partial {\overset{̅}{x}}^{2}})}^{- 1} ε^{22} & 0 \\ 0 & 0 & \frac{\partial {\tilde{x}}^{1}}{\partial {\overset{}{x}}^{1}} \frac{\partial {\tilde{x}}^{2}}{\partial {\overset{̅}{x}}^{2}} ε^{33} \end{array}) .

Owing to the increased resolution at an interface the occurring derivatives of x͂^m with respect to x̄^m are close to zero in the direction x̄^m. Formulating the application of factorization rules for Eq. (39), all terms are multiplied with such a minimizing function before the Fourier transform is carried out. For instance,

{〈 {\overset{̅}{ε}}^{11} 〉}_{21}^{ASR} = {〚 \frac{\partial {\tilde{x}}^{1}}{\partial {\overset{̅}{x}}^{1}} {({〈 {\overset{̅}{ε}}^{11} 〉}_{2}^{ASR})}^{- 1} 〛}_{1}^{- 1},

Hence, any jump discontinuities at the interfaces of the terms ε^mn are reduced, no matter if they are caused by the transformation to matched coordinates or the permittivity or permeability itself. Therefore, an increased convergence can be expected, because the over- and undershoots occurring at jump discontinuities are proportional to the jump height. An appropriate transformation for the example with the cylinders is as follows:

{\tilde{x}}^{m} ({\overset{̅}{x}}^{m}) = {\begin{array}{c} - a_{1} {({\overset{̅}{x}}^{m})}^{2} + a_{-} {\overset{̅}{x}}^{m} for 0 \leq {\overset{̅}{x}}^{m} < {\overset{̅}{x}}_{-}^{m} \\ a_{2} + a_{3} {\overset{̅}{x}}^{m} + a_{4} \sin \frac{2 π ({\overset{̅}{x}}^{m} - {\overset{̅}{x}}_{-}^{m})}{{\overset{̅}{x}}_{+}^{m} - {\overset{̅}{x}}_{+}^{m}} for {\overset{̅}{x}}_{-}^{m} \leq {\overset{̅}{x}}^{m} \leq {\overset{̅}{x}}_{+}^{m}, \\ a_{1} {({\overset{̅}{x}}^{m})}^{2} - a_{+} {\overset{̅}{x}}^{m} + a_{5} for {\overset{̅}{x}}_{+}^{m} < {\overset{̅}{x}}^{m} \leq P \end{array}

with

a_{\pm} = 2 \frac{{\tilde{x}}_{-}^{m}}{{({\overset{̅}{x}}_{-}^{m})}^{2}} {\overset{̅}{x}}_{\pm}^{m} - (1 - η) \frac{P \pm {\overset{̅}{x}}_{+}^{m}}{{\overset{̅}{x}}_{-}^{m}}, a_{1} = \frac{{\tilde{x}}_{-}^{m}}{{({\overset{̅}{x}}_{-}^{m})}^{2}} - \frac{1 - η}{{\overset{̅}{x}}_{-}^{m}},

The parameter η with 0 ≤ η ≤ 1 defines the strength of the transformation. Note that the above transformation is only valid if the interfaces are symmetric in the unit cell with P - x͂^m ₊ = x͂^m _- and P - x̄^m ₊ = x̄^m _-. The new interface positions x̄^m _± do not necessarily have to coincide with the original positions x͂^m _±.

A disadvantage of coordinate transformations in a Fourier modal method is that also uniform slabs have to be calculated numerically. To avoid problems with degenerate eigenstates we follow the procedure described in [5] and replace the solution of all far field channels by the Fourier transformation of the analytical solution in the basis of the new coordinates.

7. Numerical examples

7.1. Dielectric cylinder

A fundamental problem that provides no simple application of factorization rules is a planar periodic arrangement of cylinders. Figure 2(a) shows coordinate curves obtained by Eq. (38) and Eq. (41) to match such a cylinder in a square lattice including adaptive spatial resolution. For the sake of simplicity, μ = 1 in all further calculations. The problem of an isolated cylinder can be solved analytically by considering boundary conditions in cylindrical coordinates and results in an effective propagation constant k ^calc in the invariant direction of the cylinder [14]. If the solution of the cylinder decays exponentially outside the cylinder and the decay length is much shorter than the size of the chosen period in the Fourier modal calculation, then the interaction of different unit cells can be neglected for this mode. Hence, the propagation constant k ^calc should be a solution of the eigenvalue problem (13).

We took a structure with vacuum outside and ε = 4 inside the cylinder. The diameter is 1 μm and the unit cell size is 1.5 μm × 1.5 μm. The fundamental mode for 619.5 THz exhibits a decay of exp(-21) at a distance of 0.5 μm from the boundary of the cylinder. The calculated propagation constant is 25 1/μm. The relative deviation of the numerical solutions from the analytically derived propagation constant can be seen in Fig. 2(b) for a zigzag approximation of the cylinder by lines that are aligned vertically or horizontally in a finite uniform Cartesian grid (circles), merely matched coordinates (squares), and additionally adaptive spatial resolution (diamonds) calculated under normal incidence. The parameter for adaptive spatial resolution is η = 0.97 and the interface positions in matched coordinates remained unchanged. Because the partial derivatives of the coordinate transformation to matched coordinates are discontinuous, the results can be improved tremendously using adaptive spatial resolution. The comparison of the formulas (33) to (36) (white markers) with the formulas by Li (black markers) gives almost identical results. Both are superior compared to the simple zigzag approximation by two orders of magnitude when using adaptive spatial resolution. Therefore, we will omit the case of zigzag approximation in the following.

Fig. 2. (a) Matched coordinates including adaptive spatial resolution for a cylinder. (b) Comparison of the effective propagation constant k ³ with the calculated solution k ^calc of a dielectric cylinder. The scheme of Li is indicated by the small black markers.

Download Full Size | PDF

Another test of the scheme is the energy conservation. In the example of a dielectric cylinder the absorption should be zero. The absorption accuracy for the cylinder with a height of 50 nm can be seen in Fig. 3(a). Obviously, the energy conservation is violated for small truncation orders, but the absorption accuracy reaches already values below 10^-3 with adaptive spatial resolution for reasonable truncation orders starting from 361 harmonics.

The deviation between the application of factorization rules according to Li and the symmetry-preserving formulas (33) to (36) is almost zero (see Fig. 3(b)). In the case of absorption including adaptive spatial resolution it is smaller than 10^-5 for all truncation orders. Obviously, the global convergence behavior is given by the Gibbs phenomenon and both schemes factorize correctly.

Fig. 3. (a) Absorption accuracy of a dielectric cylinder to show energy conservation using the symmetry-preserving formulas (white markers) and Li (black markers). (b) Differences in absorption accuracy A and the effective propagation constant k ³ between the scheme by Li and our formulation.

Download Full Size | PDF

7.2. Metallic cylinder

A more complicated problem than dielectric cylinders are systems including metals. Then, a correct description of the geometry is necessary for a high accuracy. The spectra of a periodic system (square lattice with periods P ₁ = P ₂ = 700 nm) including 50 nm thin cylinders with a diameter of 300 nm consisting of a Drude metal (gold) with the plasma frequency ω _Plasma = 1.37×10¹⁶s^-1 and a damping frequency ω_γ = 0.85×10¹⁴s^-1 can be seen in Fig. 4(a) for normal incidence. Here, the same transformation as for the dielectric cylinder with adaptive spatial resolution has been used. Note that no scheme exists so far to predict the optimum choice of parameters for adaptive spatial resolution and η = 0.97 has been simply taken as a standard value of similar rectangular metallic systems. Hence, it can be that another value leads to better results, but it is not our intention to discuss this problem in closer details now.

The resonance around 370 THz has been chosen to compare the convergence of the transmission intensity with and without adaptive spatial resolution (Fig. 4(b)). The convergence can be improved for both schemes of factorization rules by adaptive spatial resolution (diamonds) compared to merely matched coordinates (squares).

The approach by Li [13] is indicated by the small black markers and exhibits a very similar behavior as the formulas (33) to (36) with the white markers (see Fig. 5(a) for details). However, the difference between two independent polarizations (P-polarization and S-polarization with incidence electric field in direction x ¹ and x ², respectively) in Fig. 5(b) shows that the result depends less on the incidence polarization when using our approach.

8. Conclusion

The presented scheme to apply factorization rules to Maxwell’s equations solves the problem of asymmetry in the original formulation [13] and exhibit the same energy conservation behavior. With the method of matched coordinates including adaptive spatial resolution it is possible to calculate complex structures efficiently in a rigorous way. The transformation to matched coordinates requires one coordinate line to be aligned along the interfaces but does not necessarily have to be continuously differentiable nor does the other coordinate line have to be normal at the interfaces. This simplifies the derivation of appropriate transformations and provides the opportunity to analyze structures that differ strongly from the unit cell. The validity of our scheme has been demonstrated by a comparison with analytical results for dielectric cylinders.

Fig. 4. (a) Spectral behavior (transmission T, reflection R, and absorption A) of an array of metallic cylinders (see the text for a detailed description of parameters) for 625 harmonics. (b) Convergence of transmission intensity at 370 THz with matched coordinates (squares) and additionally adaptive spatial resolution (diamonds).

Download Full Size | PDF

Fig. 5. (a) Difference in transmission between the scheme by Li and our formulation. (b) Spurious transmission asymmetry between two independent polarizations P and S for the totally symmetric cylinder at the resonance.

Download Full Size | PDF

Acknowledgments

We acknowledge support form BMBF (13N 9155, 13N 10146), DFG (FOR 557, FOR 730), ANR Chair of Excellence Program, the Russian Academy of Sciences, and the Russian Foundation for Basic Research.

References

1. D. M. Whittaker and I. S. Culshaw, “Scattering Matrix treatment of patterned multilayer photonic structures,” Phys. Rev. B 60, 2610–2618 (1999). [CrossRef]

2. S. G. Tikhodeev, A. L. Yablonskii, E. A. Muljarov, N. A. Gippius, and T. Ishihara, “Quasiguided modes and optical properties of photonic crystal slabs,” Phys. Rev. B 66, 45,102–1–17 (2002). [CrossRef]

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

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

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

6. P. Lalanne, “Improved formulation of the coupled-wave method for two-dimensional gratings,” J. Opt. Soc. Am. A 14, 1592–1598 (1997). [CrossRef]

7. E. Popov and M. Neviére, “Maxwell equations in Fourier space: fast-converging formulation for diffraction by arbitrary shaped, periodic, anisotropic media,” J. Opt. Soc. Am. A 18, 2886–2894 (2001). [CrossRef]

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

9. P. Götz, T. Schuster, K. Frenner, S. Rafler, and W. Osten, “Normal vector method for the RCWA with automated vector field generation,” Opt. Ex. 16, 17,295–17,301 (2008).

10. J. H. Heinbockel, “Introduction to Tensor Calculus and Continuum Mechanics,” (1996).

11. E. J. Post, Formal structure of Electromagnetism (North Holland, Amsterdam, 1962).

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

14. J. D. Jackson, Classical Electrodynamics, 3rd ed. (John Wiley & Sons, 1999).

Matched coordinates and adaptive spatial resolution in the Fourier modal method

Abstract

1. Introduction

2. Coordinate transformations

3. Covariant formulation of Maxwell’s equations

4. Factorization rules

5. Application of factorization rules to Maxwell’s equations

6. Matched coordinates and adaptive spatial resolution

7. Numerical examples

7.1. Dielectric cylinder

7.2. Metallic cylinder

8. Conclusion

Acknowledgments

References

Cited By

Figures (5)

Equations (45)

Optics Express